1. setwd and env

library(ggplot2)
library(reshape2)
#library(Rmisc)
library(Hmisc)
library(EnhancedVolcano)
library(pracma)
library(car)
library(ggrepel)
library(edgeR)
library(grid)
library(gridExtra)
library(Mfuzz)
library(M3C)
library(preprocessCore)
library(rlist)
library(ggsci)
library(scales)
library(data.table)
library(RColorBrewer)
library(readxl)
#setwd("/home/ligh/github/code_forpublication/macaca_multiple_tissue/")
source('./subroutines.R')
source('./subroutines_for_MCMT_aging.R')

2. load data

#whole body data
load('./data/pro.whole.fdr0.01_v20210108_from_NOVO_remap_solid_tissues.Rdata')
load('./data/mrna.whole_v20210108_solid_tissues.Rdata')
load('./data/met_whole_from_novo_v20210108_solid_tissues.Rdata')

# tissue data
load('./data/pro.tissues_v20210108_solid_tissues.Rdata')
load('./data/mrna.tissues_v20210108_solid_tissues.Rdata')
load('./data/met.tissues_v20210108_solid_tissues_new.Rdata')


load('./data/met.header.all.hmdb_curated.Rdata')
idx = !is.na(met.header.all.hmdb$hmdbid_highconfidence)
met.header.all.hmdb.v = met.header.all.hmdb[idx,]

3. quality control

# protein remove outliners
tmp = pro.whole
pp = prcomp(t(tmp),cor=F)
outlinerids = c()
thetissues = unique(pro.whole.info$tissue_en)
for(i in 1:length(thetissues)){
    idx = pro.whole.info$tissue_en == thetissues[i]
    tmpinfo = pro.whole.info[idx,]
    outx = is.outliner(pp$x[idx,1])
    outy = is.outliner(pp$x[idx,2])
    #outy = is.outliner(pp$x[idx,2],coef = 3)
    if(sum(outx | outy) > 0){
        outlinerids = c(outlinerids,rownames(tmpinfo)[outx | outy])
    } 
}
outlinerids
##  [1] "X06080_Skin_of_back"   "X11062_Pituitary"      "X94356_Liver"         
##  [4] "X12092_Thyroid_gland"  "X06080_Thyroid_gland"  "X94356_Cecum"         
##  [7] "X16002_Adrenal_gland"  "X92338_Adrenal_gland"  "X12390_Fallopian_tube"
## [10] "X16068_Hypothalamus"   "X06070_Pancreas"       "X16086_Uterus"
vid = !is.element(colnames(pro.whole),outlinerids)
pro.whole = pro.whole[,vid]
pro.whole.info = pro.whole.info[vid,]
# met remove outliners
tmp = met.whole
pp = prcomp(t(tmp),cor=F)
outlinerids = c()
thetissues = unique(met.whole.info$tissue_en)
for(i in 1:length(thetissues)){
    idx = met.whole.info$tissue_en == thetissues[i]
    tmpinfo = met.whole.info[idx,]
    outx = is.outliner(pp$x[idx,1])
    outy = is.outliner(pp$x[idx,2])
    #outy = is.outliner(pp$x[idx,2],coef = 3)
    if(sum(outx | outy) > 0){
        outlinerids = c(outlinerids,rownames(tmpinfo)[outx | outy])
    } 
}
outlinerids
## NULL
vid = !is.element(colnames(met.whole),outlinerids)
met.whole = met.whole[,vid]
met.whole.info = met.whole.info[vid,]
# mrna remove outliners
tmp = mrna.whole
pp = prcomp(t(tmp),cor=F)
outlinerids = c()
thetissues = unique(mrna.whole.info$tissue_en)
for(i in 1:length(thetissues)){
    idx = mrna.whole.info$tissue_en == thetissues[i]
    tmpinfo = mrna.whole.info[idx,]
    outx = is.outliner(pp$x[idx,1])
    outy = is.outliner(pp$x[idx,2])
    #outy = is.outliner(pp$x[idx,2],coef = 3)
    if(sum(outx | outy) > 0){
        outlinerids = c(outlinerids,rownames(tmpinfo)[outx | outy])
    } 
}
outlinerids
## NULL
vid = !is.element(colnames(mrna.whole),outlinerids)
mrna.whole = mrna.whole[,vid]
mrna.whole.info = mrna.whole.info[vid,]

4. set colors

alltissues = names(pro.tissues)

tissue.systems = c('Integumentary','Endocrine','Brain','Respiratory','Digestive',
             'Cardiovascular','Cardiovascular','Brain','Digestive','Endocrine',
             'Cardiovascular','Muscle','Reproductive','Digestive','Brain',
             'Immune','Renal','Endocrine','Digestive','Reproductive',
             'Digestive','Brain','Muscle','Immune','Integumentary','Endocrine',
            'Brain','Immune','Cardiovascular','Reproductive')
names(tissue.systems) = alltissues
tissue.color = pal_npg()(10)
names(tissue.color) = levels(factor(tissue.systems))

tissue.systems
##            Skin_of_back               Pituitary            Frontal_pole 
##         "Integumentary"             "Endocrine"                 "Brain" 
##                    Lung                   Liver        Arteria_cruralis 
##           "Respiratory"             "Digestive"        "Cardiovascular" 
##            Femoral_vein             Hippocampus               Ileocecum 
##        "Cardiovascular"                 "Brain"             "Digestive" 
##           Thyroid_gland         Arteria_carotis                  Muscle 
##             "Endocrine"        "Cardiovascular"                "Muscle" 
##                   Ovary                   Cecum Superior_temporal_gyrus 
##          "Reproductive"             "Digestive"                 "Brain" 
##                  Spleen                  Kidney           Adrenal_gland 
##                "Immune"                 "Renal"             "Endocrine" 
##                Duodenum          Fallopian_tube                 Stomach 
##             "Digestive"          "Reproductive"             "Digestive" 
##            Hypothalamus                   Heart                  Thymus 
##                 "Brain"                "Muscle"                "Immune" 
##             Facial_skin                Pancreas     Supramarginal_gyrus 
##         "Integumentary"             "Endocrine"                 "Brain" 
##                 Adipose             Aortic_arch                  Uterus 
##                "Immune"        "Cardiovascular"          "Reproductive"
tissue.color
##          Brain Cardiovascular      Digestive      Endocrine         Immune 
##    "#E64B35FF"    "#4DBBD5FF"    "#00A087FF"    "#3C5488FF"    "#F39B7FFF" 
##  Integumentary         Muscle          Renal   Reproductive    Respiratory 
##    "#8491B4FF"    "#91D1C2FF"    "#DC0000FF"    "#7E6148FF"    "#B09C85FF"
mypal = pal_aaas()(10)
mypal
##  [1] "#3B4992FF" "#EE0000FF" "#008B45FF" "#631879FF" "#008280FF" "#BB0021FF"
##  [7] "#5F559BFF" "#A20056FF" "#808180FF" "#1B1919FF"
show_col(mypal)

5. Figure 1: overall data distribution

5.1 Figure 1A flowchart

5.2 Figure 1B

num_omics = data.frame(num_mrna = rep(0,length(alltissues)),stringsAsFactors = F,
                       tissues = alltissues,
                       tissue_systems = tissue.systems,
                      num_protein = rep(0,length(alltissues)),
                       num_met = rep(0,length(alltissues))
                       )
for(i in 1:length(alltissues)){
    num_omics$num_mrna[i] = nrow(mrna.tissues[[alltissues[i]]])
    num_omics$num_protein[i] = nrow(pro.tissues[[alltissues[i]]])
    num_omics$num_met[i] = nrow(met.tissues[[alltissues[i]]])
}
rownames(num_omics) = alltissues
idx = sort.int(num_omics$num_protein,decreasing = F,index.return = T)$ix
num_omics.v = num_omics[idx,]
idx = sort.int(num_omics.v$tissue_systems,decreasing = F,index.return = T)$ix
num_omics.v = num_omics.v[idx,]
writetxt(num_omics.v[,c(2,3,1,4,5)],'./out/20230217_aging/Figure1_Overall_data_distribution/Figure_1B_number_omics_data.txt')
idx = sort.int(num_omics$num_protein,decreasing = F,index.return = T)$ix
num_omics.v = num_omics[idx,]
idx = sort.int(num_omics.v$tissue_systems,decreasing = F,index.return = T)$ix
num_omics.v = num_omics.v[idx,]

pomics = list()
pomics[[1]] = ggplot(num_omics.v,aes(x = factor(tissues,level = tissues),
                                     y = num_met,
                       color = tissue_systems,fill = tissue_systems))+ 
  geom_bar(stat="identity")+
theme_classic()+lghplot.addtheme(legend.position = 'none',hjust = 1,size = 10)+ 
  scale_color_aaas()+scale_fill_aaas()+ 
 theme(axis.text.x=element_blank(),axis.title.x=element_blank())+ 
  ylab('Number of metabolites')


pomics[[2]] = ggplot(num_omics.v,aes(x = factor(tissues,level = tissues),
                                     y = num_protein,
                       color = tissue_systems,fill = tissue_systems))+ 
  geom_bar(stat="identity")+
theme_classic()+lghplot.addtheme(legend.position = 'none',hjust = 1,size = 10)+
  scale_color_aaas()+scale_fill_aaas()+
ylab('Number of proteins')+ xlab('')

#pdf(file = './out/20230217_aging/Figure1_Overall_data_distribution/
#    Figure1B_number_of_omicsV__promet.pdf',width = 7,height = 5)
grid.arrange(arrangeGrob(grobs = pomics,ncol = 1,
                         top = 'Identified proteins and metabolites',
                         heights = c(3.2,5.8)))

#dev.off()
# total identified molecules
vmrna = c()
for (i in 1:length(mrna.tissues)){
    vmrna = unique(c(vmrna,rownames(mrna.tissues[[i]])))
}
length(vmrna)
## [1] 16614
vproteins = c()
for (i in 1:length(pro.tissues)){
    vproteins = unique(c(vproteins,pro.tissues.header[[i]]$Gene))
}
length(vproteins)
## [1] 14295
vmets = c()
for (i in 1:length(met.tissues)){
    vmets = unique(c(vmets,rownames(met.tissues[[i]])))
}
length(vmets)
## [1] 5789

5.3 Figure 1C_1E

5.3.1 Figure 1C mRNA

### for mRNA
mrna.whole.std = standardise_matrix(mrna.whole)

p = tsne(mrna.whole.std,labels=tissue.systems[mrna.whole.info$tissues],
         legendtextsize = 10,dotsize = 2)
p = p + theme_classic()+lghplot.addtheme(legend.position = 'none')+ 
  scale_color_aaas()+ xlab('tsne 1') + ylab('tsne 2')+
theme(axis.line = element_line(size = 1.0)) + ggtitle('Transcriptome')
#pdf(file ="./out/20230217_aging/Figure1_Overall_data_distribution/
#    Figure1C_tSNE_mrna_using_standardised.pdf",height = 4,width = 4)
print(p)

#dev.off()

5.3.2 Figure 1D protein

ida = rowSums(pro.whole < 0.1) < 0.2 * ncol(pro.whole)
sum(ida)
## [1] 3087
pro.whole.cons = pro.whole[ida,]
pro.whole.std = standardise_matrix(pro.whole.cons)

p = tsne(pro.whole.std,labels=tissue.systems[pro.whole.info$tissue_en],
         legendtextsize = 10,dotsize = 2)
p = p + theme_classic()+lghplot.addtheme(legend.position = 'none')+ 
  scale_color_aaas()+ xlab('tsne 1') + ylab('tsne 2') + ggtitle('Proteome')
#pdf(file ="./out/20230217_aging/Figure1_Overall_data_distribution/
#    Figure1D_tSNE_protein_using_standardised.pdf",height = 4,width = 4)
print(p)

#dev.off()

5.3.3 Figure 1E metabolism

met.whole.std = standardise_matrix(met.whole)
p = tsne(met.whole.std,labels=tissue.systems[met.whole.info$Tissues],
         legendtextsize = 10,dotsize = 2)
p = p + theme_classic()+lghplot.addtheme(legend.position = 'none')+ 
  scale_color_aaas()+ xlab('tsne 1') + ylab('tsne 2')+
theme(axis.line = element_line(size = 1.0)) + ggtitle('Metabolome')
#pdf(file ="./out/20230217_aging/Figure1_Overall_data_distribution/
#    Figure1E_tSNE_met_using_std.pdf",height = 4,width = 4)
print(p)

#dev.off()

6 Figure 2: tissue aging DEGs and GO

# sub routine plot volcano
plot_Volcano <- function(res2, title){
    

  tmpup = res2$Pvalue
  tmpup[is.na(tmpup)] = 1
  tmpup[res2$log2FC < 0] = 1
  sortid_up = sort.int(-log10(tmpup),decreasing = T,index.return = T)$ix
  tmpid.up = res2$ID[sortid_up[1:5]]

  tmpdown = res2$Pvalue
  tmpdown[is.na(tmpdown)] = 1
  tmpdown[res2$log2FC > 0] = 1
  sortid_down = sort.int(-log10(tmpdown),decreasing = T,index.return = T)$ix
  tmpid.down = res2$ID[sortid_down[1:5]]

  vid = is.element(res2$ID,c(tmpid.up,tmpid.down))
  tlab = res2$ID
  tlab[!vid] = NA

  keyvals <- rep('gray50', nrow(res2))
  names(keyvals) <- rep('NS', nrow(res2))   

  keyvals[which(res2$log2FC > 0.58 & res2$Pvalue < 0.05)] <- "Brown"    
  names(keyvals)[which(res2$log2FC > 0.58 & res2$Pvalue < 0.05)] <- 'High'    

  keyvals[which(res2$log2FC < -0.58 & res2$Pvalue < 0.05)] <- "darkblue"    
  names(keyvals)[which(res2$log2FC < -0.58 & res2$Pvalue < 0.05)] <- 'Low'   
  p = EnhancedVolcano(res2,
                        lab = tlab,
                        x = 'log2FC',
                        y = 'Pvalue',
                        caption = NULL,
                        title = title,
                        border = 'full',
                        titleLabSize = 12,
                        FCcutoff = 0.58,
                        cutoffLineWidth = F,
                        axisLabSize = 12,
                        subtitle = NULL,
                        cutoffLineCol = 'white',
                        gridlines.minor = F,
                        gridlines.major = F,
                        xlab = bquote(~Log[2]~ 'fold change'),
                        pCutoff = 0.05,
                        colCustom = keyvals,
                        colAlpha = 4/5,
                        legendPosition = 'none',
                        legendLabSize = 5,
                        legendIconSize = 3,
                        drawConnectors = TRUE,
                        widthConnectors = 0.5,
                        pointSize = -0.3*log10(res2$Pvalue),labSize = 3,
                        colConnectors = 'black')
    return(p)
}


list_to_matrix <- function(DEproFC,alltissues){
    DEproFC_matrix = list()
    for(i in 1:length(alltissues)){
        tmp = matrix(DEproFC[[i]],1,length(DEproFC[[i]]))
        tmp = as.data.frame(tmp)
        colnames(tmp) = names(DEproFC[[i]])
        DEproFC_matrix[[i]] = tmp
    }
    DEproFC_matrix = t(as.matrix(rbindlist(DEproFC_matrix,fill = T)))
    colnames(DEproFC_matrix) = names(DEproFC)
    #vid = rowSums(is.na(DEproFC_matrix)) < ncol(DEproFC_matrix)/2
    #DEproFC_matrix = DEproFC_matrix[vid,]
    return(DEproFC_matrix)
}

met.class_enrichment <- function(mets,annote){
  require(clusterProfiler)

  vmet = intersect(mets,rownames(annote))
  fluxgmt = data.frame(ont = annote$sub_class,
                       gene = rownames(annote),stringsAsFactors = F) 
  
  Recon3D <- enricher(gene = vmet,
                TERM2GENE=fluxgmt,
                pAdjustMethod = "BH",
                minGSSize = 1,
                pvalueCutoff =1,
                qvalueCutoff = 1
                )
  Recon3Dout = Recon3D@result
  Recon3Dout$DB = rep('HMDBclass',dim(Recon3Dout)[1])
  n = dim(Recon3Dout)[2]
  Recon3Dout = Recon3Dout[,c(n,1:(n-1))]
  return(Recon3Dout)
}

6.1 protein

DEproFC = list()
DEproPvalue = list()
DEproAging = list()
DEproFC.develop = list()
DEproPvalue.develop = list()
for(i in 1:length(alltissues)){
    #Mfuzz
    thistissue = alltissues[i]
    thispro = pro.tissues[[thistissue]]
    thispro.header = pro.tissues.header[[thistissue]]
    thispro = delete_dup_genes_forprotein(thispro,pro.tissues.header[[thistissue]])
    thispro.info = pro.tissues.info[[thistissue]]
    thisDEpro = DEGenes.simplified(thispro,catagory = thispro.info$stage == 4,
                    subset = thispro.info$stage == 4 | thispro.info$stage == 1)
    thisDEpro.develop = DEGenes.simplified(thispro,catagory = thispro.info$stage == 2,
                    subset = thispro.info$stage == 2 | thispro.info$stage == 1)
    forsort = thisDEpro$Pvalue
    forsort[is.na(forsort)] = 1
    idx = sort.int(forsort,decreasing = F,index.return = T)$ix
    DEproAging[[i]] = thisDEpro[idx,-5]
    DEproAging[[i]]$log2FC.devControl = thisDEpro.develop$log2FC
    DEproAging[[i]]$Pvalue.devControl = thisDEpro.develop$Pvalue
    
    DEproFC[[i]] = thisDEpro$log2FC
    names(DEproFC[[i]]) = rownames(thisDEpro)
    DEproPvalue[[i]] = thisDEpro$Pvalue
    names(DEproPvalue[[i]])  = rownames(thisDEpro)
    
    DEproFC.develop[[i]] = thisDEpro.develop$log2FC
    names(DEproFC.develop[[i]]) = rownames(thisDEpro.develop)
    DEproPvalue.develop[[i]] = thisDEpro.develop$Pvalue
    names(DEproPvalue.develop[[i]])  = rownames(thisDEpro.develop)
}
names(DEproAging) = alltissues
names(DEproFC) = alltissues
names(DEproPvalue) = alltissues
names(DEproFC.develop) = alltissues
names(DEproPvalue.develop) = alltissues
proteinVoconoPlot = list()
names.DEproAging = names(DEproAging)
for(i in 1:length(DEproAging)){
    res2 = DEproAging[[i]]
    
    proteinVoconoPlot[[i]] = plot_Volcano(res2,names.DEproAging[i])
}
pdf(file = './out/20230217_aging/Figure2_DEG_GO_tissue/
    Figure2A_DEpro_each_tissueV1.pdf',width = 3*5+1,height =3*6)
grid.arrange(arrangeGrob(grobs = proteinVoconoPlot,ncol = 5))
dev.off()
## png 
##   2
#Data S2
openxlsx::write.xlsx(DEproAging, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S2_DE_tissue_Aging_pro.xlsx")

# reducce size Data S2
tpath = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S2_DE_tissue_Aging_pro.xlsx"
sheetNames = openxlsx::getSheetNames(tpath)
xx = list()
for(i in 1:length(sheetNames)){
    tmp = openxlsx::read.xlsx(tpath,sheet = sheetNames[i])
    tmp[c(2,3,4,7,8)] = signif(tmp[c(2,3,4,7,8)],3)
    tmp = tmp[,-c(5,6)]
    xx[[i]] = tmp
}
names(xx) = sheetNames
openxlsx::write.xlsx(xx, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/reduce_Data S2_DE_tissue_Aging_pro.xlsx")
# num up and down proteins
DEproFC_matrix = list_to_matrix(DEproFC,alltissues)
DEproPvalue_matrix = list_to_matrix(DEproPvalue,alltissues)

DEproFC_matrix.develop = list_to_matrix(DEproFC.develop,alltissues)
DEproPvalue_matrix.develop = list_to_matrix(DEproPvalue.develop,alltissues)

Aging_pro_sigup_matrix = (DEproFC_matrix > 0.58 & DEproPvalue_matrix < 0.05  
                          & DEproPvalue_matrix.develop > 0.05) +0
Aging_pro_sigdown_matrix = -((DEproFC_matrix < -0.58 & DEproPvalue_matrix < 0.05 
                              & DEproPvalue_matrix.develop > 0.05) +0)
Aging_pro_sigall_matrix = Aging_pro_sigup_matrix + Aging_pro_sigdown_matrix

Aging_pro_updown = data.frame(stringsAsFactors = F,
                       num.up = colSums(Aging_pro_sigup_matrix,na.rm =T)/
                          colSums(!is.na(Aging_pro_sigup_matrix)),
                       num.down = colSums(Aging_pro_sigdown_matrix,na.rm =T)/
                          colSums(!is.na(Aging_pro_sigdown_matrix)),
                       num.all = colSums(abs(Aging_pro_sigall_matrix),na.rm =T)/
                         colSums(!is.na(Aging_pro_sigall_matrix)),
                        tissues = colnames(Aging_pro_sigup_matrix),
                        tissue_systems = tissue.systems)
rownames(Aging_pro_updown) = Aging_pro_updown$tissues

Aging_pro_updown
##                              num.up    num.down    num.all
## Skin_of_back            0.036514823 -0.02422270 0.06081448
## Pituitary               0.024935277 -0.02820548 0.05316973
## Frontal_pole            0.022538363 -0.01742606 0.03996803
## Lung                    0.033295982 -0.01888302 0.05220176
## Liver                   0.034220532 -0.01801802 0.05225225
## Arteria_cruralis        0.025117739 -0.04575335 0.07101526
## Femoral_vein            0.041103448 -0.03337010 0.07450331
## Hippocampus             0.013891271 -0.02023320 0.03413379
## Ileocecum               0.021394879 -0.04090909 0.06238408
## Thyroid_gland           0.037166086 -0.01646011 0.05364059
## Arteria_carotis         0.024580336 -0.02637890 0.05096942
## Muscle                  0.062205062 -0.06260720 0.12483912
## Ovary                   0.070639717 -0.08727348 0.15794957
## Cecum                   0.040690691 -0.01725431 0.05795796
## Superior_temporal_gyrus 0.015867713 -0.01687270 0.03274307
## Spleen                  0.019508449 -0.03287250 0.05238900
## Kidney                  0.018391573 -0.01220532 0.03060201
## Adrenal_gland           0.042047532 -0.01391427 0.05604055
## Duodenum                0.030563661 -0.06327753 0.09387003
## Fallopian_tube          0.029572113 -0.04311911 0.07273930
## Stomach                 0.027610674 -0.05751735 0.08519833
## Hypothalamus            0.011605416 -0.00902498 0.02063185
## Heart                   0.043579314 -0.02086957 0.06451613
## Thymus                  0.075391850 -0.32503133 0.40075259
## Facial_skin             0.009045226 -0.01404917 0.02313301
## Pancreas                0.011964948 -0.01736638 0.02933738
## Supramarginal_gyrus     0.017594835 -0.01743904 0.03504522
## Adipose                 0.008937121 -0.08154781 0.09053103
## Aortic_arch             0.066784870 -0.07934633 0.14626454
## Uterus                  0.030058007 -0.04924376 0.07935949
##                                         tissues tissue_systems
## Skin_of_back                       Skin_of_back  Integumentary
## Pituitary                             Pituitary      Endocrine
## Frontal_pole                       Frontal_pole          Brain
## Lung                                       Lung    Respiratory
## Liver                                     Liver      Digestive
## Arteria_cruralis               Arteria_cruralis Cardiovascular
## Femoral_vein                       Femoral_vein Cardiovascular
## Hippocampus                         Hippocampus          Brain
## Ileocecum                             Ileocecum      Digestive
## Thyroid_gland                     Thyroid_gland      Endocrine
## Arteria_carotis                 Arteria_carotis Cardiovascular
## Muscle                                   Muscle         Muscle
## Ovary                                     Ovary   Reproductive
## Cecum                                     Cecum      Digestive
## Superior_temporal_gyrus Superior_temporal_gyrus          Brain
## Spleen                                   Spleen         Immune
## Kidney                                   Kidney          Renal
## Adrenal_gland                     Adrenal_gland      Endocrine
## Duodenum                               Duodenum      Digestive
## Fallopian_tube                   Fallopian_tube   Reproductive
## Stomach                                 Stomach      Digestive
## Hypothalamus                       Hypothalamus          Brain
## Heart                                     Heart         Muscle
## Thymus                                   Thymus         Immune
## Facial_skin                         Facial_skin  Integumentary
## Pancreas                               Pancreas      Endocrine
## Supramarginal_gyrus         Supramarginal_gyrus          Brain
## Adipose                                 Adipose         Immune
## Aortic_arch                         Aortic_arch Cardiovascular
## Uterus                                   Uterus   Reproductive
mean(Aging_pro_updown$num.all)
## [1] 0.07529843

6.2 mRNA

DEmrnaFC = list()
DEmrnaPvalue = list()
DEmrnaAging = list()
DEmrnaFC.develop = list()
DEmrnaPvalue.develop = list()
for(i in 1:length(alltissues)){
    #Mfuzz
    thistissue = alltissues[i]
    thismrna = mrna.tissues[[thistissue]]
    #thismrna.header = mrna.tissues.header[[thistissue]]
    thismrna.info = mrna.tissues.info[[thistissue]]
    if (sum(thismrna.info$stage == 1) ==1){ 
        thismrna = cbind(thismrna[,thismrna.info$stage == 1],thismrna)
        thismrna.info = rbind(thismrna.info[thismrna.info$stage == 1,],thismrna.info)
    }
    if (sum(thismrna.info$stage == 4) ==1){ 
        thismrna = cbind(thismrna,thismrna[,thismrna.info$stage == 4])
        thismrna.info = rbind(thismrna.info,thismrna.info[thismrna.info$stage == 4,])
    }
    
    thisDEmrna = DEGenes.simplified(thismrna,catagory = thismrna.info$stage == 4,
                  subset = thismrna.info$stage == 4 | thismrna.info$stage == 1)
    
    thisDEmrna.develop = DEGenes.simplified(thismrna,catagory = thismrna.info$stage == 2,
                    subset = thismrna.info$stage == 2 | thismrna.info$stage == 1)
    
    idx = sort.int(thisDEmrna$Pvalue,decreasing = F,index.return = T)$ix
    DEmrnaAging[[i]] = thisDEmrna[idx,-5]
    DEmrnaAging[[i]]$log2FC.devControl = thisDEmrna.develop$log2FC
    DEmrnaAging[[i]]$Pvalue.devControl = thisDEmrna.develop$Pvalue
    
    
    DEmrnaFC[[i]] = thisDEmrna$log2FC
    names(DEmrnaFC[[i]]) = rownames(thisDEmrna)
    DEmrnaPvalue[[i]] = thisDEmrna$Pvalue
    names(DEmrnaPvalue[[i]])  = rownames(thisDEmrna)
    
    DEmrnaFC.develop[[i]] = thisDEmrna.develop$log2FC
    names(DEmrnaFC.develop[[i]]) = rownames(thisDEmrna.develop)
    DEmrnaPvalue.develop[[i]] = thisDEmrna.develop$Pvalue
    names(DEmrnaPvalue.develop[[i]])  = rownames(thisDEmrna.develop)
    
    
}
names(DEmrnaFC) = alltissues
names(DEmrnaPvalue) = alltissues
names(DEmrnaAging) = alltissues
names(DEmrnaFC.develop) = alltissues
names(DEmrnaPvalue.develop) = alltissues



DEmrnaFC_matrix = list_to_matrix(DEmrnaFC,alltissues)
DEmrnaPvalue_matrix = list_to_matrix(DEmrnaPvalue,alltissues)

DEmrnaFC_matrix.develop = list_to_matrix(DEmrnaFC.develop,alltissues)
DEmrnaPvalue_matrix.develop = list_to_matrix(DEmrnaPvalue.develop,alltissues)
#data S1
openxlsx::write.xlsx(DEmrnaAging, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S1_DE_tissue_Aging_mRNA.xlsx")

# reducce size Data S1
tpath = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S1_DE_tissue_Aging_mRNA.xlsx"
sheetNames = openxlsx::getSheetNames(tpath)
xx = list()
for(i in 1:length(sheetNames)){
    tmp = openxlsx::read.xlsx(tpath,sheet = sheetNames[i])
    tmp[c(2,3,4,7,8)] = signif(tmp[c(2,3,4,7,8)],3)
    tmp = tmp[,-c(5,6)]
    xx[[i]] = tmp
}
names(xx) = sheetNames
openxlsx::write.xlsx(xx, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/reduce_Data S1_DE_tissue_Aging_mRNA.xlsx")
Aging_mrna_sigup_matrix = (DEmrnaFC_matrix > 0.58 & DEmrnaPvalue_matrix < 0.05 &
                             DEmrnaPvalue_matrix.develop > 0.05) +0
Aging_mrna_sigdown_matrix = -((DEmrnaFC_matrix < -0.58 & DEmrnaPvalue_matrix < 0.05 &
                                 DEmrnaPvalue_matrix.develop > 0.05) +0)
Aging_mrna_sigall_matrix = Aging_mrna_sigup_matrix + Aging_mrna_sigdown_matrix

Aging_mrna_updown = data.frame(stringsAsFactors = F,
                               num.up = colSums(Aging_mrna_sigup_matrix,na.rm =T)/
                                 colSums(!is.na(Aging_mrna_sigup_matrix)),
                              num.down = colSums(Aging_mrna_sigdown_matrix,na.rm =T)/
                                colSums(!is.na(Aging_mrna_sigdown_matrix)),
                              num.all = colSums(abs(Aging_mrna_sigall_matrix),na.rm =T)/
                                colSums(!is.na(Aging_mrna_sigall_matrix)),
                             tissues = colnames(Aging_mrna_sigup_matrix),
                             tissue_systems = tissue.systems)
rownames(Aging_mrna_updown) = Aging_mrna_updown$tissues
Aging_mrna_updown
##                              num.up     num.down     num.all
## Skin_of_back            0.014877204 -0.013853904 0.028731108
## Pituitary               0.016756193 -0.020444159 0.037200353
## Frontal_pole            0.017589731 -0.020362887 0.037952619
## Lung                    0.008547009 -0.015562006 0.024109015
## Liver                   0.015846995 -0.011202186 0.027049180
## Arteria_cruralis        0.016141630 -0.013364575 0.029506205
## Femoral_vein            0.012257843 -0.012772158 0.025030002
## Hippocampus             0.020030204 -0.019791749 0.039821954
## Ileocecum               0.054235073 -0.028506085 0.082741158
## Thyroid_gland           0.026949541 -0.013761468 0.040711009
## Arteria_carotis         0.021807382 -0.029953331 0.051760713
## Muscle                  0.014691302 -0.028301037 0.042992339
## Ovary                   0.025969730 -0.035232818 0.061202547
## Cecum                   0.050902852 -0.040853011 0.091755862
## Superior_temporal_gyrus 0.015590200 -0.023783010 0.039373210
## Spleen                  0.029801597 -0.034164814 0.063966411
## Kidney                  0.020382901 -0.013863674 0.034246575
## Adrenal_gland           0.017783149 -0.026588398 0.044371547
## Duodenum                0.119320039 -0.126266754 0.245586793
## Fallopian_tube          0.049622774 -0.034300381 0.083923155
## Stomach                 0.009763363 -0.008439517 0.018202879
## Hypothalamus            0.001815598 -0.003473319 0.005288917
## Heart                   0.009539101 -0.011616925 0.021156026
## Thymus                  0.006846780 -0.001572909 0.008419689
## Facial_skin             0.039662514 -0.040293329 0.079955843
## Pancreas                0.024070432 -0.024445069 0.048515501
## Supramarginal_gyrus     0.014104710 -0.022870348 0.036975058
## Adipose                 0.014931669 -0.013497554 0.028429222
## Aortic_arch             0.029785810 -0.027024766 0.056810576
## Uterus                  0.020403106 -0.025028911 0.045432017
##                                         tissues tissue_systems
## Skin_of_back                       Skin_of_back  Integumentary
## Pituitary                             Pituitary      Endocrine
## Frontal_pole                       Frontal_pole          Brain
## Lung                                       Lung    Respiratory
## Liver                                     Liver      Digestive
## Arteria_cruralis               Arteria_cruralis Cardiovascular
## Femoral_vein                       Femoral_vein Cardiovascular
## Hippocampus                         Hippocampus          Brain
## Ileocecum                             Ileocecum      Digestive
## Thyroid_gland                     Thyroid_gland      Endocrine
## Arteria_carotis                 Arteria_carotis Cardiovascular
## Muscle                                   Muscle         Muscle
## Ovary                                     Ovary   Reproductive
## Cecum                                     Cecum      Digestive
## Superior_temporal_gyrus Superior_temporal_gyrus          Brain
## Spleen                                   Spleen         Immune
## Kidney                                   Kidney          Renal
## Adrenal_gland                     Adrenal_gland      Endocrine
## Duodenum                               Duodenum      Digestive
## Fallopian_tube                   Fallopian_tube   Reproductive
## Stomach                                 Stomach      Digestive
## Hypothalamus                       Hypothalamus          Brain
## Heart                                     Heart         Muscle
## Thymus                                   Thymus         Immune
## Facial_skin                         Facial_skin  Integumentary
## Pancreas                               Pancreas      Endocrine
## Supramarginal_gyrus         Supramarginal_gyrus          Brain
## Adipose                                 Adipose         Immune
## Aortic_arch                         Aortic_arch Cardiovascular
## Uterus                                   Uterus   Reproductive
mean(Aging_mrna_updown$num.all)
## [1] 0.04937392
mrnaVoconoPlot = list()
names.DEmrnaAging = names(DEmrnaAging)
for(i in 1:length(DEmrnaAging)){
    res2 = DEmrnaAging[[i]]
    
    mrnaVoconoPlot[[i]] = plot_Volcano(res2,names.DEmrnaAging[i])
}
#this plot is large, plot to file
pdf(file = './out/20230217_aging/Figure2_DEG_GO_tissue/FigureS2B_DEmrna_each_tissueV1.pdf',
    width = 3*5+1,height =3*6)
grid.arrange(arrangeGrob(grobs = mrnaVoconoPlot,ncol = 5))
dev.off()
## png 
##   2

6.3 metabolism

DEmetFC = list()
DEmetPvalue = list()
DEmetAging = list()
DEmetFC.develop = list()
DEmetPvalue.develop = list()
for(i in 1:length(alltissues)){
    #Mfuzz
    thistissue = alltissues[i]
    thismet = met.tissues[[thistissue]]
    thismet.header = met.tissues.header[[thistissue]]
    thismet.info = met.tissues.info[[thistissue]]
    if (sum(thismet.info$stage == 1) ==1){ 
        thismet = cbind(thismet[,thismet.info$stage == 1],thismet)
        thismet.info = rbind(thismet.info[thismet.info$stage == 1,],thismet.info)
    }
    if (sum(thismet.info$stage == 4) ==1){ 
        thismet = cbind(thismet,thismet[,thismet.info$stage == 4])
        thismet.info = rbind(thismet.info,thismet.info[thismet.info$stage == 4,])
    }
    
    thisDEmet = DEGenes.simplified(thismet,catagory = thismet.info$stage == 4,
                               subset = thismet.info$stage == 4 | thismet.info$stage == 1)
    
    thisDEmet.develop = DEGenes.simplified(thismet,catagory = thismet.info$stage == 2,
                               subset = thismet.info$stage == 2 | thismet.info$stage == 1)
    
    forsort = thisDEmet$Pvalue
    forsort[is.na(forsort)] = 1
    idx = sort.int(forsort,decreasing = F,index.return = T)$ix
    DEmetAging[[i]] = thisDEmet[idx,-5]
    DEmetAging[[i]]$log2FC.devControl = thisDEmet.develop$log2FC
    DEmetAging[[i]]$Pvalue.devControl = thisDEmet.develop$Pvalue
    
    DEmetFC[[i]] = thisDEmet$log2FC
    names(DEmetFC[[i]]) = rownames(thisDEmet)
    DEmetPvalue[[i]] = thisDEmet$Pvalue
    names(DEmetPvalue[[i]])  = rownames(thisDEmet)
    
    DEmetFC.develop[[i]] = thisDEmet.develop$log2FC
    names(DEmetFC.develop[[i]]) = rownames(thisDEmet.develop)
    DEmetPvalue.develop[[i]] = thisDEmet.develop$Pvalue
    names(DEmetPvalue.develop[[i]])  = rownames(thisDEmet.develop)
    
    
}
names(DEmetAging) = alltissues
names(DEmetFC) = alltissues
names(DEmetPvalue) = alltissues
DEmetFC_matrix = list_to_matrix(DEmetFC,alltissues)
DEmetPvalue_matrix = list_to_matrix(DEmetPvalue,alltissues)

DEmetFC_matrix.develop = list_to_matrix(DEmetFC.develop,alltissues)
DEmetPvalue_matrix.develop = list_to_matrix(DEmetPvalue.develop,alltissues)
#write Data S3
openxlsx::write.xlsx(DEmetAging, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S3_DE_tissue_Aging_metabolite.xlsx")

# reducce size Data S2
tpath = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S3_DE_tissue_Aging_metabolite.xlsx"
sheetNames = openxlsx::getSheetNames(tpath)
xx = list()
for(i in 1:length(sheetNames)){
    tmp = openxlsx::read.xlsx(tpath,sheet = sheetNames[i])
    tmp[c(2,3,4,7,8)] = signif(tmp[c(2,3,4,7,8)],3)
    tmp = tmp[,-c(5,6)]
    xx[[i]] = tmp
}
names(xx) = sheetNames
openxlsx::write.xlsx(xx, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/reduce_Data S3_DE_tissue_Aging_metabolite.xlsx")
# num DEmets
Aging_met_sigup_matrix = (DEmetFC_matrix > 0.58 & DEmetPvalue_matrix < 0.05 &
                            DEmetPvalue_matrix.develop > 0.05) + 0
Aging_met_sigdown_matrix = -((DEmetFC_matrix < -0.58 & DEmetPvalue_matrix < 0.05  &
                                DEmetPvalue_matrix.develop > 0.05) + 0)
Aging_met_sigall_matrix = Aging_met_sigup_matrix + Aging_met_sigdown_matrix

Aging_met_updown = data.frame(stringsAsFactors = F,
                              num.up = colSums(Aging_met_sigup_matrix,na.rm =T)/
                                colSums(!is.na(Aging_met_sigup_matrix)),
                              num.down = colSums(Aging_met_sigdown_matrix,na.rm =T)/
                                colSums(!is.na(Aging_met_sigdown_matrix)),
                              num.all = colSums(abs(Aging_met_sigall_matrix),na.rm =T)/
                                colSums(!is.na(Aging_met_sigall_matrix)),
                             tissues = colnames(Aging_met_sigup_matrix),
                             tissue_systems = tissue.systems)
rownames(Aging_met_updown) = Aging_met_updown$tissues
Aging_met_updown
##                              num.up     num.down    num.all
## Skin_of_back            0.036030341 -0.027180784 0.06321113
## Pituitary               0.127226463 -0.078032231 0.20525869
## Frontal_pole            0.013303769 -0.047302291 0.06060606
## Lung                    0.015033073 -0.021046302 0.03607937
## Liver                   0.006362059 -0.029496819 0.03585888
## Arteria_cruralis        0.020116807 -0.020765737 0.04088254
## Femoral_vein            0.000000000 -0.011150235 0.01115023
## Hippocampus             0.014792899 -0.015532544 0.03032544
## Ileocecum               0.012651265 -0.058305831 0.07095710
## Thyroid_gland           0.006051437 -0.084720121 0.09077156
## Arteria_carotis         0.048211509 -0.010886470 0.05909798
## Muscle                  0.011284722 -0.015625000 0.02690972
## Ovary                   0.040852575 -0.080521018 0.12137359
## Cecum                   0.147406266 -0.083204931 0.23061120
## Superior_temporal_gyrus 0.065759637 -0.077853364 0.14361300
## Spleen                  0.053291536 -0.035736677 0.08902821
## Kidney                  0.009730967 -0.062392673 0.07212364
## Adrenal_gland           0.021821632 -0.013757116 0.03557875
## Duodenum                0.034172662 -0.050959233 0.08513189
## Fallopian_tube          0.006607930 -0.042951542 0.04955947
## Stomach                 0.006823821 -0.047766749 0.05459057
## Hypothalamus            0.030780781 -0.024024024 0.05480480
## Heart                   0.028863500 -0.023451594 0.05231509
## Thymus                  0.023887588 -0.136299766 0.16018735
## Facial_skin             0.002056555 -0.011825193 0.01388175
## Pancreas                0.026440678 -0.016949153 0.04338983
## Supramarginal_gyrus     0.126301180 -0.045801527 0.17210271
## Adipose                 0.002421308 -0.007263923 0.00968523
## Aortic_arch             0.015799868 -0.048716261 0.06451613
## Uterus                  0.031384615 -0.035692308 0.06707692
##                                         tissues tissue_systems
## Skin_of_back                       Skin_of_back  Integumentary
## Pituitary                             Pituitary      Endocrine
## Frontal_pole                       Frontal_pole          Brain
## Lung                                       Lung    Respiratory
## Liver                                     Liver      Digestive
## Arteria_cruralis               Arteria_cruralis Cardiovascular
## Femoral_vein                       Femoral_vein Cardiovascular
## Hippocampus                         Hippocampus          Brain
## Ileocecum                             Ileocecum      Digestive
## Thyroid_gland                     Thyroid_gland      Endocrine
## Arteria_carotis                 Arteria_carotis Cardiovascular
## Muscle                                   Muscle         Muscle
## Ovary                                     Ovary   Reproductive
## Cecum                                     Cecum      Digestive
## Superior_temporal_gyrus Superior_temporal_gyrus          Brain
## Spleen                                   Spleen         Immune
## Kidney                                   Kidney          Renal
## Adrenal_gland                     Adrenal_gland      Endocrine
## Duodenum                               Duodenum      Digestive
## Fallopian_tube                   Fallopian_tube   Reproductive
## Stomach                                 Stomach      Digestive
## Hypothalamus                       Hypothalamus          Brain
## Heart                                     Heart         Muscle
## Thymus                                   Thymus         Immune
## Facial_skin                         Facial_skin  Integumentary
## Pancreas                               Pancreas      Endocrine
## Supramarginal_gyrus         Supramarginal_gyrus          Brain
## Adipose                                 Adipose         Immune
## Aortic_arch                         Aortic_arch Cardiovascular
## Uterus                                   Uterus   Reproductive
mean(Aging_met_updown$num.all)
## [1] 0.07502263
metVoconoPlot = list()
names.DEmetAging = names(DEmetAging)
for(i in 1:length(DEmetAging)){
    res2 = DEmetAging[[i]]
    
    metVoconoPlot[[i]] = plot_Volcano(res2,names.DEmetAging[i])
}

pdf(file = './out/20230217_aging/Figure2_DEG_GO_tissue/
    FigureS2C_DEmet_each_tissueV1.pdf',width = 3*5+1,height =3*6)
grid.arrange(arrangeGrob(grobs = metVoconoPlot,ncol = 5))
dev.off()
## png 
##   2

6.4 Figure 2A perc changed mols each tissue

fmt_dcimals <- function(decimals=0){
    function(x) format(x,nsmall = decimals,scientific = FALSE)
}
plotAgingNum = list()
idx = sort.int(Aging_pro_updown$num.all,decreasing = F,index.return = T)$ix
Aging_pro_updown.v = Aging_pro_updown[idx,]
idx = sort.int(Aging_pro_updown.v$tissue_systems,decreasing = F,index.return = T)$ix
Aging_pro_updown.v = Aging_pro_updown.v[idx,]
tissueindex = rownames(Aging_pro_updown.v)

#metabolites
Aging_met_updown.v = Aging_met_updown[tissueindex,]

p1 = ggplot(Aging_met_updown.v,aes(x = factor(tissues,level = tissues),y = num.up,
                       color = tissue_systems,fill = tissue_systems))+ 
  geom_bar(stat="identity",alpha = 0.8)

plotAgingNum[[1]] = p1+ geom_bar(stat="identity",aes(y = num.down),alpha = 0.6)+
  geom_hline(yintercept = 0)+
theme_classic()+lghplot.addtheme(legend.position = 'none',hjust = 1,size = 13.5)+
theme(axis.text.y = element_text(size = 9.5, face = "bold", color = "black"))+
  scale_color_aaas(alpha = 0.6)+
scale_fill_aaas(alpha = 0.6)+theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),axis.line.x =element_blank())+
scale_y_continuous(labels = fmt_dcimals(2))+
ylab('Metabolites')+ xlab('')

#protein
Aging_pro_updown.v = Aging_pro_updown[tissueindex,]

p1 = ggplot(Aging_pro_updown.v,aes(x = factor(tissues,level = tissues),y = num.up,
                       color = tissue_systems,fill = tissue_systems))+ 
  geom_bar(stat="identity",alpha = 0.8)

plotAgingNum[[2]] = p1+ geom_bar(stat="identity",aes(y = num.down),alpha = 0.6)+
  geom_hline(yintercept = 0)+
theme_classic()+lghplot.addtheme(legend.position = 'none',hjust = 1,size = 13.5)+
theme(axis.text.y = element_text(size = 9.5, face = "bold", color = "black"))+
  scale_color_aaas(alpha = 0.6)+
scale_fill_aaas(alpha = 0.6)+  
scale_y_continuous(labels = fmt_dcimals(2))+
ylab('Proteins')+ xlab('')

#mrna
Aging_mrna_updown.v = Aging_mrna_updown[tissueindex,]

p1 = ggplot(Aging_mrna_updown.v,aes(x = factor(tissues,level = tissues),y = num.up,
                       color = tissue_systems,fill = tissue_systems))+ 
  geom_bar(stat="identity",alpha = 0.8)

plotAgingNum[[3]] = p1+ geom_bar(stat="identity",aes(y = num.down),alpha = 0.6)+
  geom_hline(yintercept = 0)+
theme_classic()+lghplot.addtheme(legend.position = 'none',hjust = 1,size = 13.5)+
theme(axis.text.y = element_text(size = 9.5, face = "bold", color = "black"))+
  scale_color_aaas(alpha = 0.6)+
scale_fill_aaas(alpha = 0.6)+scale_y_continuous(labels = fmt_dcimals(2))+
ylab('mRNAs')+ xlab('')

#pdf(file = './out/20230217_aging/Figure2_DEG_GO_tissue/
#   Figure2A_number_of_Aging_moleculars_prometv1_aaas.pdf',width = 12,height =5)
grid.arrange(arrangeGrob(grobs = plotAgingNum[1:2],
                         top = 'Overall molecular changes with aging',
                         ncol = 1,heights = c(1.7,3.3)))

#dev.off()

6.5 Figure 2B GO enrichment

#write for metascape up
inputGOenrich = data.frame(tissues = names(DEproAging),stringsAsFactors = F,
                      genes = rep('c',length(DEproAging)))
#colnames(inputGOenrich) = c('tissue','GeneSymbol')
for(i in 1:length(DEproAging)){
    tmp = DEproAging[[i]]
    tgenes = unique(rownames(tmp)[tmp$Pvalue < 0.05 & tmp$log2FC > 0.58 &
                                    tmp$Pvalue.devControl > 0.05])
    tgenes = tgenes[!is.na(tgenes)]
    tgenes = tgenes[tgenes != '']
    tgenes = paste0(tgenes,collapse =',')
    inputGOenrich$genes[i] = tgenes
}
writetxt(inputGOenrich,'./out/20230217_aging/Figure2_DEG_GO_tissue/protein/_inputGOenrich_age_DEpro_up.txt',row.names = F,col.names = F)


#write for metascape down
inputGOenrich = data.frame(tissues = names(DEproAging),stringsAsFactors = F,
                      genes = rep('c',length(DEproAging)))
#colnames(inputGOenrich) = c('tissue','GeneSymbol')
for(i in 1:length(DEproAging)){
    tmp = DEproAging[[i]]
    tgenes = unique(rownames(tmp)[tmp$Pvalue < 0.05 & tmp$log2FC < -0.58 &
                                    tmp$Pvalue.devControl > 0.05])
    tgenes = tgenes[!is.na(tgenes)]
    tgenes = tgenes[tgenes != '']
    tgenes = paste0(tgenes,collapse =',')
    inputGOenrich$genes[i] = tgenes
}
writetxt(inputGOenrich,'./out/20230217_aging/Figure2_DEG_GO_tissue/protein/_inputGOenrich_age_DEpro_down.txt',row.names = F,col.names = F)
outids = c()
tpath1 = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/Enrichment_heatmap/HeatmapSelectedGO.csv'
tpath2 = './out/20230217_aging/Figure2_DEG_GO_tissue/protein//metascape_DEpro_down/Enrichment_heatmap/HeatmapSelectedGO.csv'
selectGO = readxl::read_xlsx(path = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/metascape_result.xlsx',
                             sheet = 'Enrichment')
go1 = file2frame(tpath1,sep = ',',header = T,row.names =2)
rownames(go1) = paste0(go1$GO,'-',rownames(go1))
go1term = go1$GO
tpath1a = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/Enrichment_GO/GO_AllLists.csv'
go1qval = file2frame(tpath1a,sep = ',',header = T)
rownames(go1qval) = paste0(go1qval$GO,'X_LogP_',go1qval$GeneList)
go1 = abs(as.matrix(go1[,-1]))
cname = colnames(go1)
for(i in 1:nrow(go1)){
    for(j in 1:ncol(go1)){
        tmpname = paste0(go1term[i],cname[j])
        go1[i,j] = abs(go1qval[tmpname,]$Log.q.value.)
    }
}
go1[is.na(go1)] = 0

go1for_write = cbind(data.frame(GO = rownames(go1)),-go1)
outlist = list()
outlist[[1]] = go1for_write;
outlist[[2]] = selectGO;
tmpnames = sort(unique(go1qval$GeneList))
for(i in 1:length(tmpnames)){
    outlist[[i+2]] = go1qval[go1qval$GeneList == tmpnames[[i]],]
}
names(outlist) = c('GOenrichmentALL','selectGO',tmpnames)
openxlsx::write.xlsx(outlist, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S4_GO_tissue_upregulated_protein.xlsx")


selectGO = readxl::read_xlsx(path = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_down/metascape_result.xlsx',sheet = 'Enrichment')
go2 = file2frame(tpath2,sep = ',',header = T,row.names =2)
rownames(go2) = paste0(go2$GO,'-',rownames(go2))
go2term = go2$GO
tpath2a = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_down/Enrichment_GO/GO_AllLists.csv'
go2qval = file2frame(tpath2a,sep = ',',header = T)
rownames(go2qval) = paste0(go2qval$GO,'X_LogP_',go2qval$GeneList)
go2 = abs(as.matrix(go2[,-1]))
cname = colnames(go2)
for(i in 1:nrow(go2)){
    for(j in 1:ncol(go2)){
        tmpname = paste0(go2term[i],cname[j])
        go2[i,j] = -abs(go2qval[tmpname,]$Log.q.value.)
    }
}
go2[is.na(go2)] = 0

go1for_write = cbind(data.frame(GO = rownames(go2)),go2)
outlist = list()
outlist[[1]] = go1for_write;
outlist[[2]] = selectGO;
tmpnames = sort(unique(go2qval$GeneList))
for(i in 1:length(tmpnames)){
    outlist[[i+2]] = go2qval[go2qval$GeneList == tmpnames[[i]],]
}
names(outlist) = c('GOenrichmentALL','selectGO',tmpnames)
openxlsx::write.xlsx(outlist, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S5_GO_tissue_downregulated_protein.xlsx")


thisgo = rbind(go1,go2)
thisgo.matrix = as.matrix(thisgo)
colnames(thisgo.matrix) = capitalize(gsub('X_LogP_','',colnames(thisgo.matrix)))
thisgo.matrix[abs(thisgo.matrix) < -log10(0.05)] = 0
writetxt(thisgo.matrix,'./out/20230217_aging/Figure2_DEG_GO_tissue/Figure2_heatmap_GOenrichment_metascape_DEpro_eachtissue_fdr.txt',row.names = T)
thisgo.matrix[thisgo.matrix > 4] = 4
thisgo.matrix[thisgo.matrix < -4] = -4

rownames(thisgo.matrix)[rownames(thisgo.matrix) == 'R-HSA-9716542-Signaling by Rho GTPases, Miro GTPases and RHOBTB3'] = 'R-HSA-9716542-Signaling by Rho GTPases'
rownames(thisgo.matrix)[rownames(thisgo.matrix) == 'R-HSA-1428517-The citric acid (TCA) cycle and respiratory electron transport'] = 'R-HSA-1428517-TCA cycle and respiratory electron transport'

graphics.off()
pheatmap::pheatmap(thisgo.matrix,cluster_rows = T,cluster_cols = T,
                   main ='Enrichment of proteome',
                   fontsize_row = 11,fontsize_col = 11,fontsize = 14,
                   treeheight_row = 20,treeheight_col = 20,legend = T,
                  color=colorRampPalette(c('#3B4992','gray99','#BB0021FF'))(50),
                  height = 10,width = 10
                  #file ="./out/20230217_aging/Figure2_DEG_GO_tissue/Figure2_heatmap_GOenrichment_metascape_DEpro_eachtissue.v1_aaas.pdf"
                   )

7 Figure 3 clustering aging type

7.1 mfuzz promet tissues

tissues = names(pro.tissues)
promet.tissues = list()
promet.tissues.info = list()
promet.tissues.Z = list()
mfuzz.promet.tissues = list()
promet.mstd.eset = list()
for(i in 1:length(tissues)){
#for(i in 1:1){
    tt = tissues[i]
    thispro = pro.tissues[[tt]]
    thispro = delete_dup_genes_forprotein(thispro,pro.tissues.header[[tt]])
    thispro.header = pro.tissues.header[[tt]]
    thispro = thispro[rowSums(is.na(thispro)) < 1/3*ncol(thispro), ]
    thismet = met.tissues[[tt]]
    rownames(thismet) = paste0('met_',rownames(thismet))
    thispromet = rbind2(thispro,thismet)
    thisinfo = pro.tissues.info[[tt]]
    thisinfo = thisinfo[colnames(thispromet),]
    thispromet.median = t(aggregate(t(thispromet), by=list(thisinfo$stage), 
                                    FUN=median, na.rm = T))
    mstd = standardise_matrix(thispromet.median)
    promet.tissues.Z[[i]] = mstd
    promet.tissues[[i]] = thispromet
    promet.tissues.info[[i]] = thisinfo
    mstd.v = mstd[rowSums(is.na(mstd))  == 0,]
    mstd.eset = new("ExpressionSet",exprs = mstd.v)
    promet.mstd.eset[[i]] = mstd.eset
    mfuzz.promet.tissues[[i]] = mfuzz(mstd.eset, c = 8,m = 1.5)

}
names(mfuzz.promet.tissues) = tissues
names(promet.tissues) = tissues
names(promet.tissues.Z) = tissues
names(promet.tissues.info) = tissues
names(promet.mstd.eset) = tissues
## not run

tissues = names(promet.tissues.mstd.eset)
for(i in 1:length(mfuzz.promet.tissues)){
    tpath = paste0('./out/20210428_aging/promet/tissues/',tissues[i])
    dir.create(tpath)
    pdf(paste0(tpath,'/mfuzz_plot_8A_promet.pdf'),width = 8,height = 4)
    mfuzz.plot(promet.tissues.mstd.eset[[i]],mfuzz.promet.tissues[[i]],
    mfrow=c(2,4),time.labels = c("Juvenile", "Young_adult","Middle_aged",
                                   "Elderly"),new.window = F)
    dev.off()
}
promet.tissues.mstd.eset = promet.mstd.eset
promet.tissues.Z.t = list()
for(i in 1:length(promet.tissues.Z)){
    tmp = as.data.frame(t(promet.tissues.Z[[i]]))
    promet.tissues.Z.t[[i]] = tmp
}
names(promet.tissues.Z.t) = names(promet.tissues.Z)

promet.whole.Z = t(as.matrix(as.data.frame(rbindlist(promet.tissues.Z.t,fill = T))))
colnames(promet.whole.Z) = paste0(rep(names(promet.tissues.Z),each = 4),'_',rep(1:4,times = 30))
promet.whole.Z.info = data.frame(tissue = rep(names(promet.tissues.Z),each = 4),
                                 stringsAsFactors = F,
                                stage = rep(1:4,times = 30))
rownames(promet.whole.Z.info) = colnames(promet.whole.Z)

7.2 Figure 3A whole body mfuzz

#this is used to reproduce the figures in manuscript
load('./out/promet_outdata.Rdata')
promet.whole.Z.v = promet.whole.Z[rowSums(is.na(promet.whole.Z)) < 0.5*120,]
promet.whole.Z.mean = t(aggregate(t(promet.whole.Z.v), 
                      by=list(promet.whole.Z.info$stage), FUN=mean, na.rm = T))
promet.whole.Z.mean  =promet.whole.Z.mean[-1,]
promet.whole.Z.mean.v = promet.whole.Z.mean[rowSums(is.na(promet.whole.Z.mean))  == 0,]
promet.whole.Z.mean.eset = new("ExpressionSet",exprs = promet.whole.Z.mean.v)
dim(promet.whole.Z.mean.v)
## [1] 5331    4
sum(substr(rownames(promet.whole.Z.mean.v),1,4) == 'met_')
## [1] 1221
#mfuzz.promet.whole = mfuzz(promet.norm.eset.stand, c = 8,m = 1.5)
mfuzz.plot(promet.whole.Z.mean.eset,mfuzz.promet.whole,mfrow=c(2,4),
           time.labels = c("Juvenile", "Young_adult","Middle_aged",
                                   "Elderly"),new.window = F)

7.3 Figure 3B Go enrichment

outids = c()
tpath = paste0('./out/20230217_aging/Figure3_trajactory_analysis/cluster1-8_metascape/metascape_result_curated.xlsx')
    my_data <- read_excel(tpath, sheet = "Enrichment")
    idx = regexpr('Summary',my_data$GroupID) >0
    outids = c(outids,my_data$Term[idx])
outids
##  [1] "R-HSA-2262752" "GO:0006412"    "GO:0045055"    "R-HSA-199991" 
##  [5] "GO:0006163"    "GO:0072594"    "GO:0006091"    "R-HSA-72203"  
##  [9] "WP3888"        "GO:0022411"    "GO:0007005"    "GO:0006520"   
## [13] "GO:0010256"    "GO:0097435"    "R-HSA-1280215" "GO:0051640"   
## [17] "GO:0060627"    "GO:1903827"    "GO:0006914"    "ko04144"
tpath = paste0('./out/20230217_aging/Figure3_trajactory_analysis/cluster1-8_metascape/Enrichment_GO/','GO_membership.csv')
thisgo = file2frame(tpath,sep = ',')
thisgo= thisgo[!duplicated(thisgo$GO),]
xid = is.element(thisgo$GO,outids)
thisgo = thisgo[xid,]

#rownames(thisgo) = paste0(thisgo$GO,':',thisgo$Description)
rownames(thisgo) = capitalize(paste0(thisgo$Description))
thisgo.matrix = as.matrix(thisgo[,substr(colnames(thisgo),1,6)== 'X_LogP'])
writetxt(thisgo.matrix,'./out/20230217_aging/Figure3_trajactory_analysis/Table S1_GOenrichment_metascape_promet_clusters.txt',row.names = T)

thisgo.matrix[thisgo.matrix > -3] = 0
colnames(thisgo.matrix) = capitalize(gsub('X_LogP_','',colnames(thisgo.matrix)))
colnames(thisgo.matrix) = gsub('Cluster','C',colnames(thisgo.matrix))

pheatmap::pheatmap(thisgo.matrix,cluster_rows = T,cluster_cols = T,
                   main = 'Enrichment analysis of each cluster',
                   fontsize_row = 9,fontsize_col = 10,fontsize = 10,
                   treeheight_row = 20,treeheight_col = 20,legend = T,
                  color=colorRampPalette(c('#3C5488FF','gray95','gray95'))(7),
                  breaks = c(-32,-24,-16,-8,0),
                  #file ="./out/20230217_aging/Figure3_trajactory_analysis/Figure3B_GOenrichment_metascape_prometA.pdf",
                 height = 3.5,width = 5.5
                  )

7.4 mfuzz for each tissue

7.4.1 get data

# construct data
tissue_trajectory = data.frame()
tissue_names = names(promet.tissues)
for(k in 1:8){
    tmpGeneList = names(mfuzz.promet.whole$cluster)[mfuzz.promet.whole$cluster == k]
    for(i in 1:length(promet.tissues)){
        
        M1 = promet.tissues[[i]]
        metadata.tissue = promet.tissues.info[[i]]
        idgene1 = intersect(tmpGeneList,rownames(M1))
        cc = repmat(as.matrix(rowMedians(M1[,metadata.tissue$stage == '1'],na.rm = T)),1,ncol(M1))
        tsd = repmat(as.matrix(apply(M1,1,sd,na.rm = T)),1,ncol(M1))
        
        M1.Z = (M1 - cc)/tsd
        M1.Z.v = M1.Z[idgene1,]
        M1.Z.v.mean = t(aggregate(t(M1.Z.v), 
                                  by=list(metadata.tissue$stage), FUN=median, na.rm = T))
        M1.Z.v.mean = M1.Z.v.mean[-1,]
        mean_x = colMeans(M1.Z.v.mean,na.rm  =T)
        mean_x = mean_x - mean_x[1]
        
        tmpdata2 = data.frame(expr = mean_x,stringsAsFactors = F,
                              stage = c(1,2,3,4),
                              group = factor(c('Juvenile','Young','Middle_aged','Elderly'),
                                     level = c('Juvenile','Young','Middle_aged','Elderly')),
                              cluster = rep(k,4),
                              tissue = rep(tissue_names[i],4)
                     )
        if (nrow(tissue_trajectory) <1){
            tissue_trajectory = tmpdata2
        }else{
            tissue_trajectory = rbind(tissue_trajectory,tmpdata2)
        }
    }

}
tissue_trajectory$tissue_system = tissue.systems[tissue_trajectory$tissue]
tissue_trajectory$tissue_system_color = tissue.color[tissue_trajectory$tissue_system]
# to matrix
tissue_trajectory_matrix = matrix(0,nrow(tissue_trajectory)/length(tissue_names),length(tissue_names))
for(i in 1:length(tissue_names)){
    tmptr = tissue_trajectory[tissue_trajectory$tissue == tissue_names[i],]
    tissue_trajectory_matrix[,i] = tmptr$expr
}
rownames(tissue_trajectory_matrix) = paste(tmptr$group,tmptr$cluster,sep = '_C')
colnames(tissue_trajectory_matrix) = tissue_names
tissue_tr_plot = list()
for (i in 1:8){
    idx = tissue_trajectory$cluster == i
    tissue_tr_plot [[i]] = ggplot(tissue_trajectory[idx,],
                                  aes(x= group,y = expr,group  =tissue)) +
        geom_line(aes(color = tissue_system),alpha = 0.3 ,
                  position = position_dodge(0.2),size = 2)+ scale_color_aaas()+ 
        lghplot.addtheme(hjust = 1,size = 14)+ ggtitle(paste0('Cluster ',i))+
         theme(axis.line = element_line(size = 1.2))+ xlab('')+ ylab('')+
        scale_y_continuous(labels = scales::comma_format(accuracy =0.1))
}
#pdf(file = "./out/20230217_aging/Figure3_trajactory_analysis/
#    Figure3c_trajectory_for_each_tissue_aaas.pdf",height = 9,width = 8)
grid.arrange(arrangeGrob(grobs = tissue_tr_plot,ncol = 3,heights = c(4,4,4),
            top = textGrob('Average molecular aging trajectory in 30 tissues',
                           gp=gpar(fontface="bold",  fontsize=20)),
            bottom=textGrob('Tissues', gp=gpar(fontface="bold",  fontsize=18)),
            left = textGrob('Z values', gp=gpar(fontface="bold",  fontsize=18),rot=90)))

#dev.off()

7.4.3 Figure 3D heatmap clustering

#heatmap cluster '#E64B35FF''#4DBBD5FF'
idx = substr(rownames(tissue_trajectory_matrix),1,3) != 'Juv'
dent = pheatmap::pheatmap(tissue_trajectory_matrix[idx,],scale = 'row',
                    main = 'Heatmap of molecular aging trajectory in 30 tissues',
                          height = 5,width = 6,angle_col = 45,             
                          fontsize_row = 8,fontsize_col = 8,
                          treeheight_row = 20,treeheight_col = 20,
                          #color=colorRampPalette(c('#3B4992','gray95','red'))(30),
                          #file ="./out/20230217_aging/Figure3_trajactory_analysis/Figure3D_heatmap_based_trajectoryV1.pdf",
                          color=colorRampPalette(c('#3C5488FF','gray95','#E64B35FF'))(30))

7.4.4 Figure 3E dendgrogam

dendcol = as.dendrogram(dent$tree_col)

labelColors = c('#3C5488FF','#E64B35FF')

clusMember = cutree(dent$tree_col,2)

# function to get color labels
colLab <- function(n) {
  if (is.leaf(n)) {
    a <- attributes(n)
    labCol <- labelColors[clusMember[which(names(clusMember) == a$label)]]
    attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol)
  }
  n
}
clusDendro = dendrapply(dendcol, colLab)

dendcolsplit = cut(dendcol,12)
class1 = unlist(cut(dendcol,12)$lower[[1]])
class2 = unlist(cut(dendcol,12)$lower[[2]])
tissueType = data.frame(tissue = c(tissue_names[class1],tissue_names[class2]),
                        stringsAsFactors = F,
              type = c(rep('Type II',length(class1)),rep('Type I',length(class2))))
rownames(tissueType) = tissueType$tissue
tissueType$class = tissueType$type
tissueClass  = tissueType

#pdf('./out/20230217_aging/Figure3_trajactory_analysis/Figure3E_dendrogamV1.pdf',width = 7,height = 5)
plot(clusDendro, main = "Molecular trajactory dendgrogram")

#dev.off()
#dev.off()

7.4.5 Figure 3F trajactory pca

# pca
tissue_pca = prcomp(t(tissue_trajectory_matrix),cor=F)
perc_tissue_pca = 100*summary(tissue_pca)$importance

#pdf('./out/20230217_aging/Figure3_trajactory_analysis/Figure3F_PCA.pdf',width = 7,height = 7)
TissueType = tissueType[tissue_names,]$type
ggplot(,aes(tissue_pca$x[,1],tissue_pca$x[,2],color = TissueType)) + 
   geom_point(size =5,alpha = 0.6) + 
   theme_classic() +lghplot.addtheme(legend.position = 'top')+  
   #stat_ellipse(lwd=1,level = 0.95) +
   geom_text_repel(aes(label = tissue_names),size = 5,box.padding = 0.5,face = 'bold')+
   xlab(paste("PC1(",as.character(round(perc_tissue_pca[2,1],1)),'%)',sep = '')) + 
   ylab(paste("PC2(",as.character(round(perc_tissue_pca[2,2],1)),'%)',sep = '')) + 
   scale_color_manual(values = c('#E64B35FF','#3C5488FF'))+
   theme(legend.text = element_text(size=12,face = 'bold'))+
   labs(title = "Molecular trajectory PCA")

#dev.off()

8 Figure 4 cluster by MAA

8.1 Figure 4A

promet.norm.eset.stand.matrix = as.matrix(promet.whole.Z.mean.eset)
variability = rep(0,8) 
for(j in 1:8){
    tgene = names(mfuzz.promet.whole$cluster)[mfuzz.promet.whole$cluster == j]
    tgene = intersect(tgene,rownames(promet.norm.eset.stand.matrix))
    bx = t(t(promet.norm.eset.stand.matrix[tgene,]) - mfuzz.promet.whole$centers[j,])
    variability[j] = mean(sqrt(rowSums(bx^2)/(ncol(bx)-1)))
}
tmpdata = data.frame(amplitude = as.vector(mfuzz.promet.whole$centers[,4]-
                                             mfuzz.promet.whole$centers[,1]),
                     stringsAsFactors = F,
                     variability = variability,
                     class = paste0('Cluster ',1:8))
#pdf(file ='./out/20230217_aging/Figure4_MAA_analysis/
#    FigureS4_cluster_amplitude_variability_8_promet.pdf',width = 7,height = 5)
ggplot(tmpdata,aes(variability,amplitude)) + geom_point(size = 4) + lghplot.addtheme(size = 18)+
  geom_text_repel(aes(label = class),size = 7)+theme(axis.line = element_line(size = 1.0))+
  ggtitle('Figure S4 cluster_amplitude_variability')

#dev.off()
tissues = names(promet.tissues.Z)
clusterdist = matrix(0,length(tissues),8)
clusteramplitude = matrix(0,length(tissues),8)
clusteramplitude_xx = matrix(0,length(tissues),8)
for(i in 1:length(tissues)){
    mstd =   promet.tissues.Z[[i]] 
    for(j in 1:8){
    tgene = names(mfuzz.promet.whole$cluster)[mfuzz.promet.whole$cluster == j]
    tgene = intersect(tgene,rownames(mstd))
    bx = t(t(mstd[tgene,]) - mfuzz.promet.whole$centers[j,])
    clusterdist[i,j] = mean(sqrt(rowSums(bx^2)/(ncol(bx)-1)),na.rm = T)
    if(length(tgene) < 2){
      clusteramplitude[i,j] = NA
      next;
    }
    tamp = colMeans(mstd[tgene,],na.rm = T)
    clusteramplitude[i,j] = abs(tamp[4]-tamp[1])
    clusteramplitude_xx[i,j] = tamp[4]-tamp[1]
  }
}
rownames(clusteramplitude) = tissues
rownames(clusteramplitude_xx) = tissues
rownames(clusterdist) = tissues

px = list()
for(j in 1:8){
  tmpdata = data.frame(amplitude = clusteramplitude_xx[,j],
                       stringsAsFactors = F,
                       variability = clusterdist[,j],
                       class = tissues,
                      tissue.systems = tissue.systems)
    px[[j]] = ggplot(tmpdata,aes(variability,amplitude,color = tissue.systems)) + 
      geom_point(size = 6) + lghplot.addtheme(size = 14)+
    scale_color_aaas()+scale_fill_aaas()+
    geom_text_repel(aes(label = class),size = 3,box.padding = 0.5)+
    theme(axis.line = element_line(size = 1.5))+ ggtitle(paste0('Cluster ',j))+ 
      xlab('') + ylab('')
}

#pdf(file = "./out/20230217_aging/Figure4_MAA_analysis/
#    FigureXS_cluster8_tissues_amplitude_variability_allV2.pdf",height = 9,width = 16)
grid.arrange(arrangeGrob(grobs = px,ncol = 4,margin = c(1,1,1,1),
                top = textGrob('Molecular alteration amplitude of each cluster in 30 tissues', 
                                         gp=gpar(fontface="bold",  fontsize=24)),
                         bottom=textGrob('Molecular alteration variability', 
                                         gp=gpar(fontface="bold",  fontsize=24)),
                         #bottom = 'mRNA expression(log2 FPKM)',
                        left = textGrob('Molecular alteration amplitude(MAA)', 
                                        gp=gpar(fontface="bold",  fontsize=24),rot=90)))

#dev.off()

8.2 Figure 4B

tmpaa = clusteramplitude_xx
colnames(tmpaa) = paste0("C",1:8)
dent.MAA = pheatmap::pheatmap(tmpaa,scale = 'none',height = 4,width = 3.5,
                          main  = 'Aging types in solid tissues using MAA',
                          fontsize_row = 7,fontsize_col = 7,
                          treeheight_row = 20,treeheight_col = 20,
                          color=colorRampPalette(c('#3B4992','gray95','red'))(30),
             #filename = './out/20230217_aging/Figure4_MAA_analysis             /FigureX_promet_heatmap_tissue_aging_byMAA.pdf'
             )

8.3 Figure 4C

dendrow = as.dendrogram(dent.MAA$tree_row)

labelColors = c('#3C5488FF','#E64B35FF')

clusMember = cutree(dent.MAA$tree_row,2)

# function to get color labels
colLab <- function(n) {
  if (is.leaf(n)) {
    a <- attributes(n)
    labCol <- labelColors[clusMember[which(names(clusMember) == a$label)]]
    attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol)
  }
  n
}
clusDendro = dendrapply(dendrow, colLab)

#pdf('./out/20230217_aging/Figure4_MAA_analysis/
#    FigureX_dendrogam_by_MAA.pdf',width = 7,height = 5)
plot(clusDendro, main = "Molecular trajactory dendgrogram by MAA")

#dev.off()
mean_clusteramplitude_xx = rowMeans(clusteramplitude_xx)
names(mean_clusteramplitude_xx) = rownames(clusteramplitude_xx)
tissues_pro = rownames(clusteramplitude_xx)
dendrow = as.dendrogram(dent.MAA$tree_row)
dendrowsplit = cut(dendrow,2.5)
class1 = unlist(cut(dendrow,2.5)$lower[[1]])
class2 = unlist(cut(dendrow,2.5)$lower[[2]])
tissueClass.MAA = data.frame(tissue = c(tissues_pro[class1],
                                        tissues_pro[class2]),stringsAsFactors = F,
                         class = c(rep('Type II',length(class1)),
                                   rep('Type I',length(class2))))
rownames(tissueClass.MAA)  = tissueClass.MAA$tissue
tissueClass.MAA$mean_aging_amplitude=mean_clusteramplitude_xx[rownames(tissueClass.MAA)]


tissueClass$class.MAA = tissueClass.MAA[rownames(tissueClass),]$class
tissueClass$mean_aging_amplitude = tissueClass.MAA[rownames(tissueClass),]$mean_aging_amplitude

tissueClass$type = rep('Undefine',nrow(tissueClass))
tissueClass$type[tissueClass$class == 'Type I' & tissueClass$class.MAA == 'Type I'] = 'Type I'
tissueClass$type[tissueClass$class == 'Type II' & tissueClass$class.MAA == 'Type II'] = 'Type II'
tissueClass
##                                          tissue     type   class class.MAA
## Pancreas                               Pancreas  Type II Type II   Type II
## Supramarginal_gyrus         Supramarginal_gyrus  Type II Type II   Type II
## Frontal_pole                       Frontal_pole  Type II Type II   Type II
## Liver                                     Liver  Type II Type II   Type II
## Hippocampus                         Hippocampus  Type II Type II   Type II
## Skin_of_back                       Skin_of_back  Type II Type II   Type II
## Cecum                                     Cecum  Type II Type II   Type II
## Lung                                       Lung  Type II Type II   Type II
## Superior_temporal_gyrus Superior_temporal_gyrus  Type II Type II   Type II
## Muscle                                   Muscle  Type II Type II   Type II
## Adrenal_gland                     Adrenal_gland  Type II Type II   Type II
## Hypothalamus                       Hypothalamus  Type II Type II   Type II
## Heart                                     Heart  Type II Type II   Type II
## Arteria_cruralis               Arteria_cruralis Undefine Type II    Type I
## Arteria_carotis                 Arteria_carotis Undefine Type II    Type I
## Adipose                                 Adipose   Type I  Type I    Type I
## Thyroid_gland                     Thyroid_gland   Type I  Type I    Type I
## Facial_skin                         Facial_skin Undefine  Type I   Type II
## Femoral_vein                       Femoral_vein Undefine  Type I   Type II
## Ileocecum                             Ileocecum   Type I  Type I    Type I
## Thymus                                   Thymus   Type I  Type I    Type I
## Aortic_arch                         Aortic_arch   Type I  Type I    Type I
## Spleen                                   Spleen   Type I  Type I    Type I
## Kidney                                   Kidney   Type I  Type I    Type I
## Ovary                                     Ovary   Type I  Type I    Type I
## Fallopian_tube                   Fallopian_tube   Type I  Type I    Type I
## Stomach                                 Stomach   Type I  Type I    Type I
## Duodenum                               Duodenum   Type I  Type I    Type I
## Pituitary                             Pituitary   Type I  Type I    Type I
## Uterus                                   Uterus   Type I  Type I    Type I
##                         mean_aging_amplitude
## Pancreas                        -0.075748775
## Supramarginal_gyrus              0.231195111
## Frontal_pole                     0.186849533
## Liver                            0.210951149
## Hippocampus                      0.044697652
## Skin_of_back                     0.153227025
## Cecum                            0.117593492
## Lung                             0.062388571
## Superior_temporal_gyrus          0.151674043
## Muscle                           0.127531951
## Adrenal_gland                    0.390895560
## Hypothalamus                     0.312324819
## Heart                            0.204160794
## Arteria_cruralis                -0.302318429
## Arteria_carotis                 -0.015271359
## Adipose                         -0.650050564
## Thyroid_gland                   -0.233352138
## Facial_skin                     -0.068819684
## Femoral_vein                     0.002050297
## Ileocecum                       -0.157174219
## Thymus                          -0.293446352
## Aortic_arch                     -0.341411484
## Spleen                          -0.253614804
## Kidney                          -0.119119447
## Ovary                           -0.148487802
## Fallopian_tube                  -0.242323566
## Stomach                         -0.151051049
## Duodenum                        -0.334290753
## Pituitary                       -0.074656135
## Uterus                          -0.150139727
require(scatterplot3d)
#typecolor = '#3B4992','gray99','#EE0000FF',
#pdf('./out/20230217_aging/Figure4_MAA_analysis/Figure4C_3d_scatterplotAA.pdf')
typecolor = tissueClass[rownames(clusteramplitude_xx),]$color


s3d <- scatterplot3d(clusteramplitude_xx[,8], clusteramplitude_xx[,2], clusteramplitude_xx[,4], grid = T,size = 1,  cex.lab = 2, cex.symbols = 4,  color=typecolor,
pch=19, 
scale.y=.75,

main="Aging amplitude",

xlab='C8',

ylab='C2 ',

zlab='C4')
s3d1 <- scatterplot3d(clusteramplitude_xx[,8], clusteramplitude_xx[,2], clusteramplitude_xx[,4], grid = T,size = 1,  cex.lab = 2, cex.symbols = 4,  # x y and z axis
              color=typecolor,

pch=19,       

scale.y=0.75,
scale.Z=1.25,

main="Aging amplitude",

xlab='C8',

ylab='C2 ',

zlab='C4',
                     

text(s3d$xyz.convert(clusteramplitude_xx[,c(8,2,4)]), labels = rownames(clusteramplitude_xx),
     cex= 1))
#dev.off()

8.4 Figure 4D

require(VennDiagram)
tmpnames = rownames(tissueClass)
Trajactory_Type_I = tmpnames[tissueClass$class == 'Type I']
Trajactory_Type_II = tmpnames[tissueClass$class == 'Type II']
MAA_Type_I = tmpnames[tissueClass$class.MAA == 'Type I']
MAA_Type_II = tmpnames[tissueClass$class.MAA == 'Type II']
venn.diagram(list(Trajactory_Type_I=Trajactory_Type_I,Trajactory_Type_II = Trajactory_Type_II,
                  MAA_Type_I = MAA_Type_I,MAA_Type_II=MAA_Type_II), 
             #fill=c("red","blue",),
             fill = c("cornflowerblue", "green", "yellow", "darkorchid1"),
             alpha=c(0.5,0.5,0.5,0.5), cex=2, cat.fontface=3, cat.cex = 1.3,margin = 0.1,
             filename="./out/20230217_aging/Figure4_MAA_analysis/Venn_trajactory_and_MAA.tiff")
## [1] 1

9 Figure 5 experimental validation

# Figure 5A all changed protein
tx = colSums(abs(Aging_pro_sigall_matrix),na.rm = T)
yy1 = tissueClass[names(tx),]
#pdf('./out/20230217_aging/Figure5_markers_HE/Figure_S5_allchanged_proteins.pdf')
ggplot(,aes(x = yy1$type,y = tx,color = yy1$type,fill = yy1$type)) +#geom_violin()+
     geom_boxplot(outlier.size = -1,alpha = 0.1,size = 0.8)+
     geom_jitter(size = 5,width = 0.3)+ theme_classic() + 
     ggtitle('# of changed proteins')+
     lghplot.addtheme()+ scale_color_manual(values=c('#EE0000','#3B4992','gray10'))+
     xlab('Tissue Type')+ylab('Number of changed proteins')

#dev.off()
wilcox.test(tx~yy1$type,subset= yy1$type != 'Undefine')     
## 
##  Wilcoxon rank sum exact test
## 
## data:  tx by yy1$type
## W = 143, p-value = 0.001914
## alternative hypothesis: true location shift is not equal to 0
tx = colSums(abs(Aging_pro_sigall_matrix),na.rm = T)/
  colSums(!is.na(Aging_pro_sigall_matrix))
yy1 = tissueClass[names(tx),]
#pdf('./out/20230217_aging/Figure5_markers_HE/Figure4B_allchanged_proteins_percentageV2.pdf')
ggplot(,aes(x = yy1$type,y = tx*100,color = yy1$type,fill = yy1$type)) +#geom_violin()+
     geom_boxplot(outlier.size = -1,alpha = 0.1,size = 0.8)+
     geom_jitter(size = 5,width = 0.3)+ theme_classic() +
     ggtitle('% of changed proteins')+
     lghplot.addtheme()+ scale_color_manual(values=c('#EE0000','#3B4992','gray10'))+
     xlab('Tissue Type')+ylab('Percentage of changed proteins(%)')

#dev.off()
wilcox.test(tx~yy1$type,subset= yy1$type != 'Undefine')
## 
##  Wilcoxon rank sum exact test
## 
## data:  tx by yy1$type
## W = 136, p-value = 0.007244
## alternative hypothesis: true location shift is not equal to 0

9.1 Figure 5B

p16 = file2frame('./data/P16_tissues_v20230419.txt')
utissue = unique(p16$Tissue)
tpvalues = data.frame(tissue = utissue,
                     p_Juvenile_vs_Elderly = rep(1,length(utissue)),
                     p_Young_vs_Elderly = rep(1,length(utissue)))
breaks_A = c(0.01,0.08,0.08,0.08,0.08,0.01,0.01,0.01)
breaks_B = c(3,3,3,3,2,2,0.06,3)
for(i in 1:length(utissue)){
    tmpdata = p16[p16$Tissue == utissue[i],]
    tmpdata$Area.Density.log2 = tmpdata$Area.Density*1e3#log2(tmpdata$Area.Density*1e3)
    tmpdata$class = factor(x = tmpdata$class,levels = c('Juvenile','Young','Elderly'))
    tmpp1 = t.test(Area.Density~class,data = tmpdata,subset = tmpdata$class != 'Young')$p.value
    tmpp2 = t.test(Area.Density~class,data = tmpdata,subset = tmpdata$class != 'Juvenile')$p.value
    tpvalues$p_Juvenile_vs_Elderly[i] = tmpp1
    tpvalues$p_Young_vs_Elderly[i] = tmpp2
    
    df2 = data.frame(tmean = aggregate(tmpdata$Area.Density.log2, 
                                       by=list(tmpdata$class), FUN=mean, na.rm = T)$x,
                tsd = aggregate(tmpdata$Area.Density.log2, by=list(tmpdata$class), 
                                FUN=sd, na.rm = T)$x,
                    class = factor(x= c('Juvenile','Young','Elderly'),
                                   levels = c('Juvenile','Young','Elderly'))
                 )
    thepath = paste0('./out/20230217_aging/Figure5_markers_HE/Figure5x_p16_stats', utissue[i],'.pdf')
    pdf(thepath)
    p = ggplot(df2, aes(x = class, y = tmean, fill = class)) + 
      geom_bar(stat = 'identity',color = 'black',position = position_dodge(),
               alpha = 0.5,size = 1.5)+
      geom_errorbar(aes(ymin = tmean,ymax = tmean+tsd),width = .3,size = 1.5)+
      lghplot.addthemeA(size = 28,sizex = 26,sizey = 26)+
      scale_fill_manual(values = c('#008B45FF','#3B4992FF','#EE0000FF'))+ 
      theme(axis.line = element_line(size = 1.2))+ xlab('')+ 
    ylab(bquote('Area Density ' ~ italic(x10) ^italic(3)))+
    #ylab('Log10 Area Density x 1e-7')+
    ggtitle(utissue[i])
    #if(i <=8){
    #    p = p+ scale_y_break(c(breaks_A[i],breaks_B[i]), scales = "free",space  = 0.2)
    #}
    
    print(p)
    dev.off()
}

tpvalues
##           tissue p_Juvenile_vs_Elderly p_Young_vs_Elderly
## 1    Aortic_arch          4.546295e-03       4.553667e-03
## 2         Spleen          2.928558e-03       2.941005e-03
## 3         Kidney          3.607333e-04       3.618602e-04
## 4          Ovary          2.399612e-03       2.487708e-03
## 5         Thymus          1.122249e-02       1.151782e-02
## 6         Uterus          1.052903e-04       1.057245e-04
## 7  Thyroid_gland          6.633705e-04       6.751286e-04
## 8        Stomach          3.601685e-05       3.616963e-05
## 9       Pancreas          2.616949e-02       5.453556e-03
## 10  Skin_of_back          2.446174e-02       4.146347e-02
## 11          Lung          4.895734e-03       6.619731e-03
## 12         Liver          2.051557e-01       2.593903e-01
## 13        Muscle          3.401860e-01       6.380951e-01

9.2 Figure 5C

p21 = file2frame('./data/P21_tissues_v20230419.txt')
utissue = unique(p21$Tissue)
tpvalues = data.frame(tissue = utissue,
                     p_Juvenile_vs_Elderly = rep(1,length(utissue)),
                     p_Young_vs_Elderly = rep(1,length(utissue)))
breaks_A = c(0.05,0.4,0.4,0.4,0.1,0.4,0.4,0.05)
breaks_B = c(1,4,4,4,2,4,4,1)
for(i in 1:length(utissue)){
    tmpdata = p21[p21$Tissue == utissue[i],]
    tmpdata$Area.Density.log2 = tmpdata$Area.Density*1e3#log2(tmpdata$Area.Density*1e3)
    tmpdata$class = factor(x = tmpdata$class,levels = c('Juvenile','Young','Elderly'))
    tmpp1 = t.test(Area.Density~class,data = tmpdata,subset = tmpdata$class != 'Young')$p.value
    tmpp2 = t.test(Area.Density~class,data = tmpdata,subset = tmpdata$class != 'Juvenile')$p.value
    tpvalues$p_Juvenile_vs_Elderly[i] = tmpp1
    tpvalues$p_Young_vs_Elderly[i] = tmpp2
    
    df2 = data.frame(tmean = aggregate(tmpdata$Area.Density.log2, 
                                       by=list(tmpdata$class), FUN=mean, na.rm = T)$x,
                tsd = aggregate(tmpdata$Area.Density.log2, by=list(tmpdata$class), 
                                FUN=sd, na.rm = T)$x,
                    class = factor(x= c('Juvenile','Young','Elderly'),
                                   levels = c('Juvenile','Young','Elderly'))
                 )
    thepath = paste0('./out/20230217_aging/Figure5_markers_HE/Figure4x_p21_stats', utissue[i],'.pdf')
    pdf(thepath)
    p = ggplot(df2, aes(x = class, y = tmean, fill = class)) + 
      geom_bar(stat = 'identity',color = 'black',position = position_dodge(),
               alpha = 0.5,size = 1.5)+
      geom_errorbar(aes(ymin = tmean,ymax = tmean+tsd),width = .3,size = 1.5)+
      lghplot.addthemeA(size = 28,sizex = 26,sizey = 26)+
      scale_fill_manual(values = c('#008B45FF','#3B4992FF','#EE0000FF'))+ 
      theme(axis.line = element_line(size = 1.2))+ xlab('')+ 
    ylab(bquote('Area Density ' ~ italic(x10) ^italic(3)))+
    #ylab('Log10 Area Density x 1e-7')+
    ggtitle(utissue[i])
    #if(i <=8){
    #    p = p+ scale_y_break(c(breaks_A[i],breaks_B[i]), scales = "free",space  = 0.2)
    #}
    
    print(p)
    dev.off()
}

tpvalues
##           tissue p_Juvenile_vs_Elderly p_Young_vs_Elderly
## 1    Aortic_arch          3.371849e-03       3.490217e-03
## 2         Spleen          6.244567e-04       6.370290e-04
## 3         Kidney          1.052679e-02       1.070538e-02
## 4          Ovary          3.484938e-04       3.559713e-04
## 5         Thymus          4.256548e-02       4.462434e-02
## 6         Uterus          2.571205e-02       2.611076e-02
## 7  Thyroid_gland          1.553183e-02       1.686152e-02
## 8        Stomach          4.794586e-05       4.974984e-05
## 9       Pancreas          2.700608e-01       1.668549e-01
## 10  Skin_of_back          5.553932e-02       2.923295e-01
## 11          Lung          6.940087e-04       6.950478e-04
## 12         Liver          9.815520e-01       9.972549e-01
## 13        Muscle          3.620859e-01       7.890578e-01

9.3 Figure 5D

Cellcounts = file2frame('./data/cell_counts_v20230407.txt')

utissue = unique(Cellcounts$Tissue)
tpvalues = data.frame(tissue = utissue,
                     p_Juvenile_vs_Elderly = rep(1,length(utissue)),
                     p_Young_vs_Elderly = rep(1,length(utissue)))
for(i in 1:length(utissue)){
    tmpdata = Cellcounts[Cellcounts$Tissue == utissue[i],]
    tmpdata$class = factor(x = tmpdata$class,
                           levels = c('Juvenile','Young','Elderly'))
    tmpp1 = t.test(number~class,data = tmpdata,subset = tmpdata$class != 'Young')$p.value
    tmpp2 = t.test(number~class,data = tmpdata,subset = tmpdata$class != 'Juvenile')$p.value
    tpvalues$p_Juvenile_vs_Elderly[i] = tmpp1
    tpvalues$p_Young_vs_Elderly[i] = tmpp2
    
    df2 = data.frame(tmean = aggregate(tmpdata$number, by=list(tmpdata$class), 
                                       FUN=mean, na.rm = T)$x,
                tsd = aggregate(tmpdata$number, by=list(tmpdata$class), FUN=sd, na.rm = T)$x,
                    class = factor(x= c('Juvenile','Young','Elderly'),
                                   levels = c('Juvenile','Young','Elderly'))
                 )
    thepath = paste0('./out/20230217_aging/Figure5_markers_HE/Figure4x_HE_cellcounts_', utissue[i],'_v1.pdf')
    pdf(thepath)
    p = ggplot(df2, aes(x = class, y = tmean, fill = class)) + 
      geom_bar(stat = 'identity',color = 'black',position = position_dodge(),alpha = 0.5,size = 1.5)+
      geom_errorbar(aes(ymin = tmean,ymax = tmean+tsd),width = .3,size = 1.5)+
      lghplot.addthemeA(size = 28,sizex = 26,sizey = 26)+
      scale_fill_manual(values = c('#008B45FF','#3B4992FF','#EE0000FF'))+ 
      theme(axis.line = element_line(size = 1.2))+ xlab('')+ 
      ylab('Number of Parenchymal Cells')+ggtitle(utissue[i])
    print(p)
    dev.off()
}

tpvalues
##          tissue p_Juvenile_vs_Elderly p_Young_vs_Elderly
## 1        Thymus          3.842297e-07       2.479669e-04
## 2       Stomach          4.735346e-08       4.458945e-06
## 3   Aortic_arch          8.038428e-07       2.041040e-07
## 4         Ovary          6.576220e-09       6.727973e-05
## 5        Spleen          1.724836e-06       4.997698e-01
## 6 Thyroid_gland          3.814005e-06       9.571183e-02
## 7        Kidney          2.734793e-01       3.090672e-01

10 Figure 6: Translation efficiency

10.1 ratio data construction

tissues = names(promet.tissues.Z)
clusterdist = matrix(0,length(tissues),8)
clusteramplitude = matrix(0,length(tissues),8)
clusteramplitude_xx = matrix(0,length(tissues),8)
for(i in 1:length(tissues)){
    mstd =   promet.tissues.Z[[i]] 
    for(j in 1:8){
    tgene = names(mfuzz.promet.whole$cluster)[mfuzz.promet.whole$cluster == j]
    tgene = intersect(tgene,rownames(mstd))
    bx = t(t(mstd[tgene,]) - mfuzz.promet.whole$centers[j,])
    clusterdist[i,j] = mean(sqrt(rowSums(bx^2)/(ncol(bx)-1)),na.rm = T)
    if(length(tgene) < 2){
      clusteramplitude[i,j] = NA
      next;
    }
    tamp = colMeans(mstd[tgene,],na.rm = T)
    clusteramplitude[i,j] = abs(tamp[4]-tamp[1])
    clusteramplitude_xx[i,j] = tamp[4]-tamp[1]
  }
}
rownames(clusteramplitude) = tissues
rownames(clusteramplitude_xx) = tissues
rownames(clusterdist) = tissues
ratio.tissues = list()
ratio.tissues.info = list()
pro.tissues.forRatio = list()
mrna.tissues.forRatio = list()
tissues = names(pro.tissues)
overlaptissues = intersect(names(pro.tissues),names(mrna.tissues))
for( i in 1:length(overlaptissues)){
    this_tissue = overlaptissues[i]
    thispro = delete_dup_genes_forprotein(pro.tissues[[this_tissue]],
                                          pro.tissues.header[[this_tissue]])
    thismrna = mrna.tissues[[this_tissue]]
    vcol = intersect(colnames(thispro),colnames(thismrna))
    vrow = intersect(rownames(thispro),rownames(thismrna))
    thispro = thispro[vrow,vcol]
    thismrna  = thismrna[vrow,vcol]
    pro.tissues.forRatio[[i]] = thispro
    mrna.tissues.forRatio[[i]] = thismrna
    ratio.tissues[[i]] = thispro - thismrna
    ratio.tissues.info[[i]] = pro.tissues.info[[i]][vcol,]
}
names(ratio.tissues) = overlaptissues
names(pro.tissues.forRatio) = overlaptissues
names(mrna.tissues.forRatio) = overlaptissues
names(ratio.tissues.info) = overlaptissues
ratio_amplitude = list()
pro_amplitude = list()
mrna_amplitude = list()
n = length(ratio.tissues)
ratio_out = data.frame(meanChangeRatio = rep(0,n),stringsAsFactors = F,
                fc_up_down = rep(0,n),
                tissues = names(ratio.tissues),
                meanChangeproZ = rep(0,n),
                fc_up_down_pro = rep(0,n),
                meanChangemrnaZ = rep(0,n),
                fc_up_down_mrna = rep(0,n))
for (i in 1:length(ratio.tissues)){
    bx = ratio.tissues[[i]]
    bx.info = ratio.tissues.info[[i]]
    id = bx.info$stage < 5
    bx = bx[,id]
    bx.info = bx.info[id,]
    cx = t(aggregate(t(bx), by=list(bx.info$stage), FUN=mean, na.rm = T))
    cx = as.matrix(standardise_1(new("ExpressionSet",exprs = cx)))
    cx = cx[-1,]
    ee = cx[,4]-cx[,1]
    ratio_amplitude[[i]] = ee;
    #hist(ee,30)
    ratio_out$meanChangeRatio[i] = mean(ee,na.rm = T)
    ratio_out$fc_up_down[i] = log2(sum(ee > 0,na.rm = T)/sum(ee < -0,na.rm = T))
    
    bx = pro.tissues.forRatio[[i]]
    bx.info = ratio.tissues.info[[i]]
    id = bx.info$stage < 5
    bx = bx[,id]
    bx.info = bx.info[id,]
    cx = t(aggregate(t(bx), by=list(bx.info$stage), FUN=mean, na.rm = T))
    cx = as.matrix(standardise_1(new("ExpressionSet",exprs = cx)))
    cx = cx[-1,]
    ee = cx[,4]-cx[,1]
    pro_amplitude[[i]] = ee
    #hist(ee,30)
    ratio_out$meanChangeproZ[i] = mean(ee,na.rm = T)
    ratio_out$fc_up_down_pro[i] = log2(sum(ee > 0,na.rm = T)/sum(ee < -0,na.rm = T))
    
    bx = mrna.tissues.forRatio[[i]]
    bx.info = ratio.tissues.info[[i]]
    id = bx.info$stage < 5
    bx = bx[,id]
    bx.info = bx.info[id,]
    cx = t(aggregate(t(bx), by=list(bx.info$stage), FUN=mean, na.rm = T))
    cx = as.matrix(standardise_1(new("ExpressionSet",exprs = cx)))
    cx = cx[-1,]
    ee = cx[,4]-cx[,1]
    mrna_amplitude[[i]] = ee
    #hist(ee,30)
    ratio_out$meanChangemrnaZ[i] = mean(ee,na.rm = T)
    ratio_out$fc_up_down_mrna[i] = log2(sum(ee > 0,na.rm = T)/sum(ee < -0,na.rm = T))
}
rownames(ratio_out) = names(ratio.tissues)
names(ratio_amplitude) = names(ratio.tissues)
names(pro_amplitude) = names(ratio.tissues)
names(mrna_amplitude) = names(ratio.tissues)

10.2 Figure 6A design

10.3 Figure S9 plot mRNA and protein

mRNA.mean = list()
pro.mean = list()
for(i in 1:length(pro.tissues.forRatio)){
    mRNA.mean[[i]] = rowMeans(mrna.tissues.forRatio[[i]],na.rm = T)
    pro.mean[[i]] = rowMeans(pro.tissues.forRatio[[i]],na.rm = T)
}
names(mRNA.mean) = names(mrna.tissues.forRatio)
names(pro.mean) = names(pro.tissues.forRatio)
RNA.v = list_to_matrix(mRNA.mean,names(mRNA.mean))
pro.v = list_to_matrix(pro.mean,names(pro.mean))

pp = list()
for(i in 1:ncol(pro.v)){
    idaa = !is.na(RNA.v[,i]) & !is.na(pro.v[,i])
    tcor = cor.test(RNA.v[idaa,i],pro.v[idaa,i])$estimate
    tmpdata = data.frame(xx = RNA.v[idaa,i],
                        yy = pro.v[idaa,i])
    pp[[i]] = ggplot(tmpdata,aes(x= xx,y = yy))+geom_point(size =1,) + theme_bw()+
        theme(plot.margin = margin(0.1,0.1,0.1,0.1,"cm"))+
        lghplot.addthemeA(size = 16,sizex = 16,sizey = 16)+        
        annotate(geom="text", x=4, y=25, 
                label = paste('R=',signif(tcor,3)),
               color="darkblue",size = 8,face = "italic")+
        xlab('')+ylab('')+ggtitle(colnames(pro.v)[i])
}

png(file = "./out/20230217_aging/Figure6_ratio_mrna_pro/Figure S9_mrna_vs_pro.png",
    width = 1250,height = 1500)
grid.arrange(arrangeGrob(grobs = pp, ncol = 5,
                         bottom=textGrob('mRNA expression(log2 CPM)', 
                                         gp=gpar(fontface="bold",  fontsize=22)),
                        left = textGrob('Protein abundance(log2 Peak Area)',
                                  gp=gpar(fontface="bold",  fontsize=22),rot=90)))
dev.off()
## png 
##   2

10.4 Figure S10 predict protein vs measured protein

ratio = pro.v - RNA.v
ratio = apply(ratio,1,median,na.rm =T)

# remove ratio = 0 and NA
#id = is.na(ratio) | ratio == 0
#ratio = ratio[!id]

pro.prediction = RNA.v + ratio

pp = list()
for(i in 1:ncol(pro.v)){
    idaa = !is.na(pro.prediction[,i]) & !is.na(pro.v[,i])
    tcor = cor.test(pro.prediction[idaa,i],pro.v[idaa,i])$estimate
    tmpdata = data.frame(xx = pro.prediction[idaa,i],
                        yy = pro.v[idaa,i])
    pp[[i]] = ggplot(tmpdata,aes(x= xx,y = yy))+geom_point(size =1,) + theme_bw()+
        theme(plot.margin = margin(0.1,0.1,0.1,0.1,"cm"))+
        lghplot.addthemeA(size = 16,sizex = 16,sizey = 16)+        
        annotate(geom="text", x=15, y=25, 
                label = paste('R=',signif(tcor,3)),
               color="darkblue",size = 8,face = "italic")+
        xlab('')+ylab('')+ggtitle(colnames(pro.v)[i])
}

png(file = "./out/20230217_aging/Figure6_ratio_mrna_pro/Figure S10_predPro_vs_pro.png",
    width = 1250,height = 1500)
grid.arrange(arrangeGrob(grobs = pp, ncol = 5,
                         bottom=textGrob('Predicted protein abundance(log2 Peak Area)', 
                                         gp=gpar(fontface="bold",  fontsize=22)),
                        left = textGrob('Protein abundance(log2 Peak Area)', 
                                        gp=gpar(fontface="bold",  fontsize=22),rot=90)))
dev.off()
## png 
##   2

10.5 Figure 6B density plot

systems.Color = substr(pal_aaas()(10),1,7)
names(systems.Color) = unique(tissue.systems)
systems.Color.fortissue = systems.Color[tissue.systems]
names(systems.Color.fortissue) = names(ratio_amplitude)
tissueClass$color = tissueClass$type
tissueClass$color[tissueClass$color == 'Type I'] = '#BB0021'
tissueClass$color[tissueClass$color == 'Type II'] = '#3B4992'
tissueClass$color[tissueClass$color == 'Undefine'] = '#B09C85'


# hist plot
p1 = list()
xplot <- function(ratio_amplitude,xtissue,tissueClass){
    tcolor = tissueClass[xtissue,]$color
    tmp =  ggplot(,aes(x = ratio_amplitude[[xtissue]]))+theme_classic()+
         #geom_histogram(binwidth=0.1,aes(y=..density..), colour="black", fill="white")+
         geom_vline(xintercept=0, linetype="dashed", color = "blue", size=1)+
         geom_density(alpha=.6, fill=tcolor) +lghplot.addtheme()+
      theme(axis.line = element_line(size = 1.2))+ 
      xlab('')+ylab('')+ggtitle(xtissue)
    return(tmp)
}

index = sort.int(ratio_out$meanChangeRatio,decreasing = F,index.return = T)$ix
vclass = tissueClass[rownames(ratio_out)[index],]
idx = sort.int(vclass$class,decreasing = F,index.return = T)$ix
vclass = vclass[idx,]
for(i in 1:nrow(vclass)){
    xtissue = vclass$tissue[i]
    p1[[i]] = xplot(ratio_amplitude,xtissue,tissueClass)
}
#pdf(file = './out/20230217_aging/Figure6_ratio_mrna_pro/Figure6B_Changed_ratio_distribute_based_on_promet_sortby_v2.pdf',width = 25,height = 12.5)
grid.arrange(arrangeGrob(grobs = p1,ncol = 8,
           top=textGrob('Distribution of changed translation efficiency to degradation (cTED) with aging', 
                                        gp=gpar(fontface="bold",  fontsize=30)),
            bottom=textGrob('Changed normalized TED', 
                                        gp=gpar(fontface="bold",  fontsize=30)),
            left = textGrob('Density', 
                                        gp=gpar(fontface="bold",  fontsize=30),rot=90)))

#graphics.off()

10.6 Figure 6C ratio plot

ratio_updown = zeros(length(ratio_amplitude),2)
tissues =names(ratio_amplitude)
rownames(ratio_updown) = tissues
colnames(ratio_updown) = c('pec_up','pec_down')
for(i in 1:length(tissues)){
    tmp = as.vector(ratio_amplitude[[i]])
    tmpN = sum(!is.na(tmp))
    ratio_updown[i,1] = sum(tmp > 0,na.rm =T)/tmpN
    ratio_updown[i,2] = sum(tmp < 0,na.rm =T)/tmpN
}
index = sort.int(ratio_out$meanChangeRatio,decreasing = F,index.return = T)$ix
ratio_updown = ratio_updown[index,]

tmp = melt(ratio_updown,stringsAsFactor =F)
colnames(tmp) = c('tissues','Direction','perctage')
tmp$tissues = factor(tmp$tissues,levels = rownames(vclass))
#pdf('./out/20230217_aging/Figure6_ratio_mrna_pro/
#   Figure6C_tissue_trans_effectiveness_amplitude_ratio_V2.pdf',width = 8,height = 7)
ggplot(data=tmp, aes(x=tissues, y=perctage, fill=Direction))+theme_classic()+
       geom_col(position = "fill",alpha = 0.8)+
      ylab('Relative ratio')+ xlab('')+scale_fill_manual(values=c('#BB0021', '#3B4992'))+
      lghplot.addtheme(legend.position = 'top',hjust = 1,size = 14)+
      ggtitle('Relative ratio of up- and down- TEDs')+
      theme(legend.text=element_text(size=20))

#graphics.off()

10.7 Figure 6D correlation

mean_aging_amplitude = rowMeans(clusteramplitude_xx[,])
mean_aging_amplitude = mean_aging_amplitude[rownames(ratio_out)]
tmpdata = data.frame(ratioChange = ratio_out$meanChangeRatio,stringsAsFactors = F,
                     ratioFC = ratio_out$fc_up_down,
                     amplitude = mean_aging_amplitude,
                    tissues = rownames(ratio_out))
rownames(tmpdata) = tmpdata$tissues
tmpdata$tissueclass = tissueClass[rownames(tmpdata),]$type
tmp1 = tmpdata
cor.test(tmp1$ratioChange,tmp1$amplitude)
## 
##  Pearson's product-moment correlation
## 
## data:  tmp1$ratioChange and tmp1$amplitude
## t = 7.9553, df = 28, p-value = 1.154e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6747763 0.9176362
## sample estimates:
##       cor 
## 0.8326318
#pdf('./out/20230217_aging/Figure6_ratio_mrna_pro/Figure6D_overall_translation_effective_score_vs_aging_amplitute_promet_V2.pdf',width = 7,height = 5)
p = ggplot(tmp1,aes(ratioChange,amplitude,color = tissueclass)) + 
  theme_classic()+
      geom_point(size = 6,aes(color = tissueclass)) + 
      lghplot.addtheme()+geom_smooth(color = 4,size = 0.5,method = 'lm')+
      geom_text_repel(aes(label = tissues),size = 4,box.padding = 0.5)+
      #scale_color_aaas(alpha = 0.6)+scale_fill_aaas(alpha = 0.6)+
      scale_color_manual(values=c('#BB0021', '#3B4992','#B09C85'))+
      theme(axis.line = element_line(size = 1.2))+
      xlab('Averaged changed translation effective score')+
      ylab('Aging amplitude')+
      ggtitle('Correlation between MAA and cTED')
print(p)

#dev.off()

11 Figure 7

11.1 GO data construct

#outids = c('R-HSA-72766','GO:0006413','hsa03010','ko04142','hsa04120','ko03050')
outids = c('R-HSA-72766','hsa03010','hsa04142','hsa04120','hsa03050')
#tpath = paste0('./out/20210428_aging/promet/tissues/metascape/Aging_up_genes_metascape/Enrichment_GO/','GO_membership.csv')
tpath = paste0('./out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/Enrichment_GO/','GO_membership.csv')
thisgo = file2frame(tpath,sep = ',')
thisgo= thisgo[!duplicated(thisgo$GO),]
xid = is.element(thisgo$GO,outids)
thisgo = thisgo[xid,]
#rownames(thisgo) = paste0(thisgo$GO,':',thisgo$Description)
rownames(thisgo) = paste0(thisgo$GO)
thisgo.matrix = -as.matrix(thisgo[,substr(colnames(thisgo),1,6)== 'X_LogP'])

# replace qval
tpath1a = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/Enrichment_GO/GO_AllLists.csv'
go1qval = file2frame(tpath1a,sep = ',',header = T)
rownames(go1qval) = paste0(go1qval$GO,'X_LogP_',go1qval$GeneList)
go1term = rownames(thisgo.matrix)
cname = colnames(thisgo.matrix)
for(i in 1:nrow(thisgo.matrix)){
    for(j in 1:ncol(thisgo.matrix)){
        tmpname = paste0(go1term[i],cname[j])
        thisgo.matrix[i,j] = abs(go1qval[tmpname,]$Log.q.value.)
    }
}
thisgo.matrix[is.na(thisgo.matrix)] = 0


thisgo.matrix[abs(thisgo.matrix) < 3] = 0
colnames(thisgo.matrix) = capitalize(gsub('X_LogP_','',colnames(thisgo.matrix)))
colnames(thisgo.matrix) = gsub('Cluster','C',colnames(thisgo.matrix))
thisgo.matrix.up = thisgo.matrix
tmpname = rownames(thisgo.matrix.up)
#add hsa03050 
thisgo.matrix.up = rbind(thisgo.matrix.up,rep(0,30))
rownames(thisgo.matrix.up) = c(tmpname,'hsa03050')
thisgo.matrix.up = thisgo.matrix.up[outids,]
thisgo.matrix.up
##             Adipose Adrenal_gland Aortic_arch Arteria_carotis Arteria_cruralis
## R-HSA-72766       0             0           0               0                0
## hsa03010          0             0           0               0                0
## hsa04142          0             0           0               0                0
## hsa04120          0             0           0               0                0
## hsa03050          0             0           0               0                0
##             Cecum Duodenum Facial_skin Fallopian_tube Femoral_vein Frontal_pole
## R-HSA-72766     0        0           0              0          4.3            0
## hsa03010        0        0           0              0          0.0            0
## hsa04142        0        0           0              0          0.0            0
## hsa04120        0        0           0              0          0.0            0
## hsa03050        0        0           0              0          0.0            0
##             Heart Hippocampus Hypothalamus Ileocecum Kidney Liver Lung Muscle
## R-HSA-72766   0.0           0            0         0      0   0.0    0      0
## hsa03010      0.0           0            0         0      0   0.0    0      0
## hsa04142      3.5           0            3         0      0   3.2    0      0
## hsa04120      0.0           0            0         0      0   0.0    0      0
## hsa03050      0.0           0            0         0      0   0.0    0      0
##             Ovary Pancreas Pituitary Skin_of_back Spleen Stomach
## R-HSA-72766     0        0         0            0      0       0
## hsa03010        0        0         0            0      0       0
## hsa04142        0        0         0            0      0       0
## hsa04120        0        0         0            0      0       0
## hsa03050        0        0         0            0      0       0
##             Superior_temporal_gyrus Supramarginal_gyrus Thymus Thyroid_gland
## R-HSA-72766                       0                   0      0             0
## hsa03010                          0                   0      0             0
## hsa04142                          0                   0      0             0
## hsa04120                          0                   0      0             0
## hsa03050                          0                   0      0             0
##             Uterus
## R-HSA-72766      0
## hsa03010         0
## hsa04142         0
## hsa04120         0
## hsa03050         0
outids = c('R-HSA-72766','hsa03010','hsa04142','hsa04120','hsa03050')
outids_des = c('Translation','Ribosome','Lysosome',
                               'Ubiquitin mediated proteolysis', 'Proteasome') 
tpath = paste0('./out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_down/Enrichment_GO/','GO_membership.csv')
thisgo = file2frame(tpath,sep = ',')
thisgo= thisgo[!duplicated(thisgo$GO),]
xid = is.element(thisgo$GO,outids)
thisgo = thisgo[xid,]
rownames(thisgo) = paste0(thisgo$GO)

# replace qval
tpath1a = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_down/Enrichment_GO/GO_AllLists.csv'
go1qval = file2frame(tpath1a,sep = ',',header = T)
rownames(go1qval) = paste0(go1qval$GO,'X_LogP_',go1qval$GeneList)
go1term = rownames(thisgo.matrix)
cname = colnames(thisgo.matrix)
for(i in 1:nrow(thisgo.matrix)){
    for(j in 1:ncol(thisgo.matrix)){
        tmpname = paste0(go1term[i],cname[j])
        thisgo.matrix[i,j] = abs(go1qval[tmpname,]$Log.q.value.)
    }
}
thisgo.matrix[is.na(thisgo.matrix)] = 0


thisgo.matrix = as.matrix(thisgo[,substr(colnames(thisgo),1,6)== 'X_LogP'])
thisgo.matrix[abs(thisgo.matrix) < 3 ] = 0
colnames(thisgo.matrix) = capitalize(gsub('X_LogP_','',colnames(thisgo.matrix)))
colnames(thisgo.matrix) = gsub('Cluster','C',colnames(thisgo.matrix))
rownames(thisgo.matrix)
## [1] "R-HSA-72766" "hsa03010"    "hsa03050"    "hsa04142"    "hsa04120"
thisgo.matrix.down = thisgo.matrix
thisgo.matrix.down = thisgo.matrix.down[outids,]
thisgo.matrix.down
##             Adipose Adrenal_gland Aortic_arch Arteria_carotis Arteria_cruralis
## R-HSA-72766    -5.2             0        -5.9               0             -4.3
## hsa03010        0.0             0         0.0               0              0.0
## hsa04142        0.0             0         0.0               0              0.0
## hsa04120        0.0             0         0.0               0              0.0
## hsa03050        0.0             0         0.0               0              0.0
##             Cecum Duodenum Facial_skin Fallopian_tube Femoral_vein Frontal_pole
## R-HSA-72766     0      0.0           0           -6.9            0            0
## hsa03010        0      0.0           0           -5.0            0            0
## hsa04142        0      0.0           0            0.0            0            0
## hsa04120        0      0.0           0            0.0            0            0
## hsa03050        0     -4.5           0            0.0            0            0
##             Heart Hippocampus Hypothalamus Ileocecum Kidney Liver Lung Muscle
## R-HSA-72766     0           0            0      -7.2      0    -6    0   -4.6
## hsa03010        0           0            0      -3.9      0     0    0    0.0
## hsa04142        0           0            0      -3.6      0     0    0    0.0
## hsa04120        0           0            0       0.0      0     0    0    0.0
## hsa03050        0           0            0       0.0      0     0    0    0.0
##             Ovary Pancreas Pituitary Skin_of_back Spleen Stomach
## R-HSA-72766  -7.7        0       -11         -3.8      0   -12.0
## hsa03010      0.0        0       -10          0.0      0    -4.7
## hsa04142      0.0        0         0          0.0      0    -5.7
## hsa04120      0.0        0         0          0.0      0     0.0
## hsa03050      0.0        0         0         -3.1      0    -6.7
##             Superior_temporal_gyrus Supramarginal_gyrus Thymus Thyroid_gland
## R-HSA-72766                       0                   0  -34.0             0
## hsa03010                          0                   0   -7.3             0
## hsa04142                          0                   0    0.0             0
## hsa04120                          0                   0   -7.8             0
## hsa03050                          0                   0  -25.0             0
##             Uterus
## R-HSA-72766     -5
## hsa03010         0
## hsa04142         0
## hsa04120         0
## hsa03050         0

11.2 Figure 7A

outids = c('R-HSA-72766','hsa03010','hsa04142','hsa04120','hsa03050')
outids_des = c('Translation','Ribosome','Lysosome',
                               'Ubiquitin mediated proteolysis', 'Proteasome')
thisgo.matrix.all = thisgo.matrix.up+ thisgo.matrix.down

tmpnames = rownames(tissueClass)
tmpnames = c(tmpnames[tissueClass$type == 'Type II'],tmpnames[tissueClass$type == 'Undefine'],
             tmpnames[tissueClass$type == 'Type I'])

thisgo.matrix.all = thisgo.matrix.all[,tmpnames]
thisgo.matrix.all[thisgo.matrix.all >= 3] = 1
thisgo.matrix.all[thisgo.matrix.all <= -3] = -1

tclass = data.frame(class = tissueClass[tmpnames,]$type,row.names = tmpnames)
colnames(tclass) <- c("Tissue_type")
ann_colors = list(
   Tissue_type= c('#FD60A7','gray85','#4fffF7')
    
)
names(ann_colors$Tissue_type) = c('Type I','Undefine', 'Type II')
#sort tissues
tmpname  = colnames(thisgo.matrix.all)
rownames(thisgo.matrix.all) = outids_des
pheatmap::pheatmap(thisgo.matrix.all,cluster_rows = F,annotation_colors = ann_colors,
                   annotation_col = tclass, 
                 main = 'Enrichment in mRNA translation and protein degradation related pathways',
                   cluster_cols = F,fontsize_row = 9,fontsize_col = 10,
                 fontsize = 10,treeheight_row = 20,treeheight_col = 20,legend = T,
                  color=colorRampPalette(c('#3B4992','gray99','#EE0000FF'))(6),
                  #file ="./out/20230217_aging/Figure6_heatmap_TED_machanism/Figure 6A_enrichment_metascape_promet_all_tissue_TEDs_V6.pdf",
                 height = 3,width = 10
                  )

idx = tclass$Tissue_type != 'Undefine'
tmpaa = thisgo.matrix.all[,idx]
tmpclass = tclass$Tissue_type[idx]
rownames(tmpaa)
## [1] "Translation"                    "Ribosome"                      
## [3] "Lysosome"                       "Ubiquitin mediated proteolysis"
## [5] "Proteasome"
cor.test(abs(tmpaa[1,]),(tmpclass == 'Type I')+0,method = 'spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  abs(tmpaa[1, ]) and (tmpclass == "Type I") + 0
## S = 1571, p-value = 0.01725
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.46291
cor.test(abs(tmpaa[2,]),(tmpclass == 'Type I')+0,method = 'spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  abs(tmpaa[2, ]) and (tmpclass == "Type I") + 0
## S = 1497.7, p-value = 0.01144
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.48795
cor.test(abs(tmpaa[3,]),(tmpclass == 'Type I')+0,method = 'spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  abs(tmpaa[3, ]) and (tmpclass == "Type I") + 0
## S = 3210.5, p-value = 0.6353
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.09759001
cor.test(abs(tmpaa[4,]),(tmpclass == 'Type I')+0,method = 'spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  abs(tmpaa[4, ]) and (tmpclass == "Type I") + 0
## S = 2340, p-value = 0.3273
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho 
## 0.2
cor.test(abs(tmpaa[5,]),(tmpclass == 'Type I')+0,method = 'spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  abs(tmpaa[5, ]) and (tmpclass == "Type I") + 0
## S = 2301.4, p-value = 0.2957
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2132007

11.3 Figure 7B

require(data.table)
amplitude_matrix = list()

for(i in 1:length(pro_amplitude)){
    tmp = matrix(pro_amplitude[[i]],1,length(pro_amplitude[[i]]))
    tmp = as.data.frame(tmp)
    colnames(tmp) = names(pro_amplitude[[i]])
    amplitude_matrix[[i]] = tmp
}
#amplitude_matrix = rbindlist(amplitude_matrix,fill = T)
amplitude_matrix = t(as.matrix(rbindlist(amplitude_matrix,fill = T)))
colnames(amplitude_matrix) = names(pro_amplitude)

vid = rowSums(is.na(amplitude_matrix)) < ncol(amplitude_matrix)/2
amplitude_matrix.v = amplitude_matrix[vid,]
amplitude_matrix.v[is.na(amplitude_matrix.v)] = NA
dim(amplitude_matrix.v)
## [1] 3902   30
tcor = cor(t(amplitude_matrix.v),ratio_out$meanChangeRatio,use = "pairwise")
names(tcor) = rownames(amplitude_matrix.v)
gsea_input = sort(tcor,decreasing = T)

calcp <- function(x,y){
    return(cor.test(x,y,use = "pairwise")$p.value)
}
tpvalue = apply(amplitude_matrix.v,1,calcp,ratio_out$meanChangeRatio)
tpvalue = tpvalue[names(gsea_input)]
xnames = names(tpvalue)
trans_initialA  = xnames[substr(xnames,1,3) == 'EIF']
trans_initialA = sort(trans_initialA)
trans_initialA
##  [1] "EIF1AX"   "EIF2A"    "EIF2AK2"  "EIF2B1"   "EIF2B2"   "EIF2B3"  
##  [7] "EIF2B4"   "EIF2B5"   "EIF2D"    "EIF2S1"   "EIF2S2"   "EIF3A"   
## [13] "EIF3B"    "EIF3D"    "EIF3E"    "EIF3F"    "EIF3G"    "EIF3H"   
## [19] "EIF3I"    "EIF3J"    "EIF3K"    "EIF3L"    "EIF3M"    "EIF4A1"  
## [25] "EIF4A2"   "EIF4A3"   "EIF4B"    "EIF4E"    "EIF4E2"   "EIF4EBP1"
## [31] "EIF4G1"   "EIF4G2"   "EIF4G3"   "EIF4H"    "EIF5A"    "EIF5A2"  
## [37] "EIF5B"    "EIF6"
trans_initialB  = xnames[substr(xnames,1,3) == 'EEF']
trans_initialB = sort(trans_initialB)
trans_initialB
## [1] "EEF1A1"    "EEF1A2"    "EEF1AKMT1" "EEF1B2"    "EEF1E1"    "EEF1G"    
## [7] "EEF2"
trans_initialA = c(trans_initialA,trans_initialB)
for(i in 1:length(pro_amplitude)){
    tx = pro_amplitude[[i]][trans_initialA]
    #tx= tx[!is.na(tx)]
    tmp = data.frame(genes = trans_initialA,stringsAsFactors = F,
                    amplitude = as.vector(tx))
    if(i ==1){
        tmpout = tmp;
    }else{
        tmpout = cbind(tmpout,tmp[,'amplitude'])
    }
    
}
tmpout = t(tmpout[,-1])

rownames(tmpout) = names(pro_amplitude)
colnames(tmpout) = trans_initialA
idx = colSums(is.na(tmpout)) < 3
sum(idx)
## [1] 31
#tmpout = tmpout[rowSums(is.na(tmpout)) < 3,]
tmpout = tmpout[,idx]


metadata <- tissueClass[rownames(tmpout),c('type','type')]
tclass = data.frame(class = metadata$type,row.names = rownames(metadata))
colnames(tclass) <- c("Tissue_type")

ann_colors = list(
   Tissue_type= c('#FD60A7','gray85','#4fffF7')
    
)
names(ann_colors$Tissue_type) = c('Type I','Undefine', 'Type II')

pheatmap::pheatmap(t(tmpout),annotation_col = tclass,cluster_rows = F,
                   annotation_colors = ann_colors,
                   main = 'Change in initiation and elongation factors',
                   color=colorRampPalette(c('#3B4992','gray99','#EE0000'))(30),
                    #file ="./out/20230217_aging/Figure6_heatmap_TED_machanism/Figure6B_heatmap_translation_elongationV3.pdf",
                  height = 8,width = 8
                  )  

11.4 Figure 7C

xnames = names(tpvalue)
ribosome = file2frame('./data/ribosome_proteins_from_kegg.txt',header = F)
ribosome = ribosome$V1
trans_initialA = intersect(ribosome,xnames)
trans_initialA
##  [1] "RPS2"   "RPS3"   "RPS3A"  "RPS4X"  "RPS5"   "RPS7"   "RPS9"   "RPS13" 
##  [9] "RPS14"  "RPS15"  "RPS16"  "RPS20"  "RPS21"  "RPS23"  "RPS24"  "RPS27L"
## [17] "RPS27A" "RPS29"  "FAU"    "RPSA"   "RPL4"   "RPL5"   "RPL6"   "RPL7"  
## [25] "RPL7A"  "RPL8"   "RPL12"  "RPL13"  "RPL13A" "RPL14"  "RPL18"  "RPL22" 
## [33] "RPL23"  "RPL23A" "RPL27"  "RPL28"  "RPL30"  "RPL32"  "RPL35"  "RPL38" 
## [41] "RPLP1"  "RPLP2"
for(i in 1:length(pro_amplitude)){
    tx = pro_amplitude[[i]][trans_initialA]
    #tx= tx[!is.na(tx)]
    tmp = data.frame(genes = trans_initialA,stringsAsFactors = F,
                    amplitude = as.vector(tx))
    if(i ==1){
        tmpout = tmp;
    }else{
        tmpout = cbind(tmpout,tmp[,'amplitude'])
    }
    
}
tmpout = t(tmpout[,-1])

rownames(tmpout) = names(pro_amplitude)
colnames(tmpout) = trans_initialA
idx = colSums(is.na(tmpout)) < 3
sum(idx)
## [1] 34
#tmpout = tmpout[rowSums(is.na(tmpout)) < 5,]
tmpout = tmpout[,idx]
tmpout = tmpout[rownames(tissueClass),]

metadata <- tissueClass#[rownames(tmpout),]
#metadata <- tissueClass[rownames(tmpout),c('class','class')]
tclass = data.frame(class = metadata$type,row.names = rownames(metadata))
colnames(tclass) <- c("Tissue_type")

ann_colors = list(
   Tissue_type= c('#FD60A7','gray85','#4fffF7')
    
)
names(ann_colors$Tissue_type) = c('Type I','Undefine', 'Type II')

pheatmap::pheatmap(t(tmpout),annotation_col = tclass,cluster_rows = F,
                   annotation_colors = ann_colors,
                   main = 'Change in ribosomal proteins',
                   cluster_cols = T,color=colorRampPalette(c('#3B4992','gray99','#EE0000FF'))(30),
                    #file ="./out/20230217_aging/Figure6_heatmap_TED_machanism/Figure6C_heatmap_ribosomeV3.pdf",
                  height = 8,width = 8
                  )